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Abstract 

The  enemy’s  synergistic  deployment  of  low  observable  (LO)  airborne  threats  and  stand¬ 
off/escort  jammers  may  defeat  USAF  fighter/interceptor  air-to-air  radars.  A  potential  solution 
to  this  problem  is  to  employ  an  adaptive  beamformer.  An  adaptive  beamformer  can  be  used  to 
cancel  interferences  while  allowing  desired  signals  to  be  passed  from  other  directions.  The  jamming 
signal  can  be  eliminated  by  either  1)  placing  a  null  in  the  direction  of  the  jammer,  or  2)  coherently 
cancelling  the  jamming  signal. 

Fighter  aircraft  equipped  with  X-Band  radars  having  an  adaptive  beamforming  capability 
exhibit  degraded  performance  due  to  radome  multipath  reflections  (RMR).  The  reflections  from 
the  radome  into  the  antenna  array  degrade  the  pattern  by  causing  the  formation  of  large  sidelobes. 
This  phenomenon  makes  detection  of  LO  targets  difficult  when  jamming  enters  the  high  RMR 
sidelobes.  The  combination  of  broadband  jammers  and  the  radar  operating  in  a  high-resolution 
mode  (wide  frequency  bandwidth)  requires  a  broadband  beamformer  to  solve  the  problem.  Future 
low  radar  cross  section  (RCS)  radomes  will  degrade  the  beampattern  even  more,  making  the  prob¬ 
lem  worse.  Actual  radomes  are  curved  and  cause  reflected  wavefronts  to  be  nonplanar.  Analyzing 
beamformer  performance  under  these  conditions  is  difficult,  since  the  beamformer  response  is  de¬ 
fined  as  the  amplitude  and  phase  change  experienced  by  a  complex  plane  wave  as  a  function  of 
frequency  and  location.  To  work  with  planar  wavefronts  and  for  the  sake  of  simplicity,  we  model 
the  radome  as  two  flat  plates  parallel  to  array  broadside  (one  at  each  end  of  the  array).  This 
radome  model  demonstrates  degraded  beamformer  performance  by  causing  large  sidelobes  to  form. 
This  approximates  the  characteristics  of  the  multipath  environment  caused  by  a  real  radome.  The 
ability  of  a  linearly  constrained,  minimum  variance  beamformer  to  cancel  interference  arriving  in 
the  sidelobes  is  evaluated.  The  cancellation  performance  depends  on  how  many  taps  are  employed 
and  the  source  frequency  bandwidth. 
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ANALYZING  ADAPTIVE  BEAMFORMER  PERFORMANCE 
IN  A  STABLE  MULTIPATH  ENVIRONMENT 


I.  Introduction 

1.1  Background 

The  enemy’s  synergistic  deployment  of  low  observable  (LO)  airborne  threats  and  stand¬ 
off/escort  jammers  may  defeat  USAF  fighter/interceptor  air-to-air  radars.  A  potential  solution 
to  this  problem  is  to  employ  an  adaptive  beamformer,  which  consists  of  an  antenna  array  with 
weights.  An  adaptive  beamformer  can  be  used  to  cancel  interferences  while  passing  desired  signals 
from  other  directions.  The  jamming  signal  can  be  eliminated  by  either  1)  placing  a  null  in  the 
direction  of  the  jammer,  or  2)  coherently  cancelling  the  jamming  signal. 

1.2  Problem 

The  presence  of  the  aircraft  radome  in  front  of  radar  antennas  degrades  adaptive  beamformer 
performance.  The  reflections  from  the  radome  into  the  antenna  array  degrade  the  pattern  by 
causing  the  formation  of  large  sidelobes.  This  phenomenon  makes  detection  of  LO  targets  difficult 
when  jamming  enters  the  high  RMR  lobes.  These  RMR  sidelobes  form  because  energy  reflects 
into  the  main  and  other  high  gain  beams.  If  a  jammer  is  located  in  the  direction  of  a  high  RMR 
sidelobe,  a  large  amount  of  the  beamformer  output  power  will  be  due  to  energy  reflected  into  the 
main  beam.  In  this  case  a  null  formed  in  the  direction  of  the  jammer  will  not  attenuate  interferer 
energy  which  enters  the  main  beam.  One  solution  is  to  put  a  null  in  the  main  beam.  This  is  not  a 
good  idea  because  this  causes  desired  signals  to  also  be  attenuated.  Thus,  jammer  signals  must  be 
coherently  cancelled  so  that  desired  signals  can  be  detected. 


1.3  Assumptions 


All  data  are  assumed  to  be  wide  sense  stationary  (WSS),  which  means  the  data  covariance 
matrix  is  independent  of  time.  The  data  received  by  the  array  are  assumed  to  be  zero  mean. 
All  desired  signals  are  assumed  to  be  uncorrelated  with  all  interfering  signals.  The  signals  and 
interferes  are  uncorrelated  with  the  noise,  which  is  white.  The  sources  are  assumed  to  originate 
as  points  in  the  far  field. 

The  radome  is  assumed  to  be  perfectly  smooth;  no  random  phase  is  associated  with  reflecting 
off  the  plates.  The  signal  is  only  attenuated  in  amplitude.  Only  a  single  bounce  off  the  plate  into 
the  array  is  analyzed.  All  discussion  of  frequency  is  in  terms  of  normalized  frequency. 

1.4  Scope 

The  scope  of  this  thesis  is  limited  to  modeling  the  radome  as  two  flat  plates.  The  direct  path 
signal  through  the  first  plate  and  the  first  bounce  Bignal  off  the  second  plate  are  generated  and 
used  as  input  to  a  fully  adaptive  beamformer.  The  output  power  as  a  function  of  the  number  of 
taps  is  calculated  to  determine  the  amount  of  interference  cancellation  obtained. 

1.5  Approach 

The  radome  is  modeled  to  simulate  its  effect  on  beamformer  performance.  Since  beamformer 
response  is  defined  in  terms  of  the  amplitude  and  phase  change  experienced  by  a  complex  plane  wave 
as  a  function  of  frequency  and  location,  a  radome  model  which  supports  plane  waves  is  developed. 
The  radome  is  modeled  as  two  flat  plates  with  reflection  coefficients  which  are  functions  of  angle. 
This  radome  model  is  used  to  simulate  the  RMR  environment. 

The  beamformer  to  be  analyzed  consists  of  an  equally  spaced,  linear  array  of  antenna  el¬ 
ements  along  with  weights.  The  number  of  antenna  elements  and  taps  (in  the  broadband  case) 
is  variable.  Taps  are  added  behind  each  sensor  in  the  broadband  case  to  provide  for  frequency 
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domain  filtering  over  the  signal  bandwidth.  For  a  detailed  discussion  on  nonzero  bandwidth  signals 
and  analysis  of  tapped  delay-line  arrays,  see  Compton’s  book  [9].  The  criteria  for  weight  selec¬ 
tion  is  linearly  constrained,  power  minimization,  also  referred  to  as  linearly  constrained,  minimum 
variance  (LCMV)  beamforming.  Before  evaluating  beamformers  in  a  multipath  environment,  nar¬ 
rowband  and  broadband  beamformers  are  designed  and  validated  by  comparing  results  to  those  in 
the  literature. 

Using  the  radome  model  described  above,  a  RMR  environment  is  generated  .  A  jammer 
is  placed  in  a  large  RMR  sidelobe,  so  that  energy  reflects  into  the  main  beam.  The  ability  of 
the  beamformer  to  cancel  this  jamming  signal  is  evaluated  by  observing  the  output  power  of  the 
beamformer  in  the  direction  of  the  jammer. 
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II.  Literature  Review 


This  section  presents  a  summary  of  research  in  the  area  of  adaptive  beamforming.  It  starts 
out  with  general  beamforming  concepts  and  then  reviews  data  generation  and  other  more  specific 
topics  related  to  the  subject  of  this  thesis. 

Van  Veen  and  Buckley  wrote  an  article  on  the  general  concepts  involved  in  beamforming  [45]. 
This  article  presents  the  concepts  on  linearly  constrained  minimum  variance  (LCMV)  beamform¬ 
ing  used  in  this  thesis.  In  particular,  we  use  the  generalized  sidelobe  canceller  (GSC)  form  of  the 
LCMV  beamformer.  The  GSC  is  the  most  general  form  of  adaptive  beamformer  and  has  demon¬ 
strated  better  results  than  any  of  the  specific  beamformers  such  as  the  multiple  sidelobe  canceller 
(MSC)  [46].  The  MSC  can  be  represented  as  a  specific  example  of  a  GSC.  The  LCMV  beamformer 
approach  is  preferred  since  it  doesn’t  assume  a  partially  adaptive  solution  like  the  MSC  does.  By 
starting  with  the  more  general  GSC,  a  partially  adaptive  solution  based  on  the  expected  interfer¬ 
ence  environment  can  be  determined.  Van  Veen  and  Buckley  [45]  cover  adaptive  algorithms  and 
implementation  issues  which  are  not  discussed  as  part  of  this  thesis.  Van  Veen  also  wrote  a  chapter 
in  a  book  on  minimum  variance  beamforming  [40].  This  material  is  updated  from  his  previous 
articles.  It  presents  some  new  information  which  makes  the  derivations  easier  to  follow. 

Most  of  the  research  in  the  area  of  beamforming  assumes  narrowband  signals  and  doesn’t  dis¬ 
cuss  the  details  of  broadband  beamforming.  Buckley  wrote  an  article  on  broadband  beamforming 
which  develops  models  for  broadband  sources  [5].  This  article,  along  with  those  of  Van  Veen,  pro¬ 
vides  all  the  details  necessary  to  design  narrowband  and  broadband  beamformers.  It  also  describes 
models  for  uncorrelated  broadband  signals  with  spatial  and/or  frequency  extent. 

Compton’s  article  and  book  provide  an  explanation  of  broadband  interference  cancellation 
performance  using  tapped  delay  lines  [8]  [9].  Monzingo  and  Miller  also  discuss  broadband  con¬ 
cepts  [30].  They  discuss  practical  aspects  of  adaptive  arrays  and  include  many  of  the  mathematical 
derivations. 
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Dudgeon’s  article  covers  the  fundamentals  of  digital  array  processing  [10];  it  applies  basic 
digital  signal  processing  principles  to  beamforming.  His  book  [11]  contains  multidimensional  signal 
processing  ideas  applied  to  beamforming.  It  describes  the  delay-and-sum  beamformer,  which  is  used 
in  this  thesis.  The  section  on  wavenumber-frequency  space  filtering  was  applied  to  take  advantage  of 
the  fast  Fourier  transform  (FFT)  in  computing  the  beamformer  response.  Another  book,  by  Widrow 
and  Stearns  [51],  covers  the  fundamentals  of  adaptive  signal  processing.  It  covers  basic  adaptive 
systems  and  their  applications,  including  beamforming.  Alexander’s  book  [1]  covers  adaptive  signal 
processing  in  more  detail  than  Widrow  and  Stearns,  but  doesn’t  include  beamforming  applications. 

Haykin  wrote  a  chapter  in  a  book  [24]  on  narrowband  signal  and  noise  models  for  radar  array 
processing  applications.  He  describes  how  off-diagonal  elements  in  the  spatial  correlation  matrix 
of  the  signal-in-space  vector  are  all  zero  and  the  spatial  correlation  matrix  of  the  received  signal 
vector  is  a  Toeplitz  matrix  when  the  sources  are  jointly  uncorrelated  or  noncoherent.  The  spatial 
correlation  matrix  of  the  received  signal  vector  is  non-Toeplitz  when  the  sources  are  correlated  or 
coherent.  One  method  of  decorrelating  a  desired  signal  from  an  interferer  is  based  on  restoring 
the  Toeplitz  structure  of  the  received  signal  vector  [17].  He  describes  a  narrowband  multipath 
model  in  detail.  Haykin  also  edited  a  book  [18]  which  contains  several  articles  on  array  processing 
applications  to  radar.  Nitzberg’s  recently  published  book  [32]  covers  array  signal  processing  for 
radar.  It  emphasizes  multiple  sidelobe  cancellers,  but  discusses  topics  of  equal  concern  to  linearly 
constrained  beamformers.  A  book  by  Hudson  [22]  has  several  sections  dealing  with  coherent  signals, 
including  multipath,  and  their  effects  on  cancellation  performance.  It  contains  discussions  not 
found  in  the  other  sources.  It  also  contains  some  relevant  and  useful  vector,  matrix,  and  statistics 
properties  in  the  appendices. 

Scharf’s  book  [35]  contains  compact,  narrowband  models  for  linear  combinations  of  vectors, 
including  one  for  the  data  matrix  associated  with  multisensor  arrays.  The  sections  on  unitary 
decomposition,  singular  value  decomposition,  and  low  rank  approximation  are  referred  to  in  this 
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thesis.  Shanmugan’s  book  [37]  was  used  as  a  reference  for  random  signal  principles.  His  definition 
of  correlation  coefficient,  along  with  other  random  signal  principles,  is  applied  in  this  study. 

An  article  by  Buckley  and  Griffiths  [6]  about  an  adaptive  GSC  with  derivative  constraints  was 
useful.  It  pointed  out  that  when  using  derivative  constraints  the  quiescent  beampattern  depends 
on  the  location  of  the  phase  center  of  the  array.  This  information  helped  to  solve  a  problem 
encountered  when  applying  eigenvector  constraints  to  an  adaptive  GSC.  Er  and  Cantoni  [13]  also 
wrote  an  article  on  derivative  constraints  for  broadband,  element  space,  array  processors.  In  all 
the  arrays  they  analyzed,  the  phase  center  of  the  array  was  chosen  as  the  zero  phase  reference. 

Frost  [15]  first  wrote  about  the  use  of  linear  constraints  for  the  time  domain,  element  space 
processors,  on  which  the  beamformers  in  this  thesis  are  based.  He  includes  a  lot  of  the  mathematics 
missing  from  the  other  articles.  He  also  provides  a  detailed  geometrical  interpretation  of  the 
linearly  constrained  minimum  variance  (LCMV)  problem.  Er  and  Cantoni  [14]  discuss  a  broadband 
processor  without  pre-steering  time  delays,  which  are  required  by  the  Frost  processor.  Griffiths  and 
Jim  first  described  the  GSC  approach  to  LCMV  beamforming  in  their  article  [21].  Jablon  [23]  also 
discussed  the  GSC  and  evaluated  its  steady  state  cancellation  performance. 

Adaptive  arrays  typically  have  main  lobe  and  sidelobe  shapes  differing  from  those  designed 
from  deterministic,  nonadaptive  arrays.  Three  articles  present  procedures  for  approximating  a 
desired  quiescent  beampattern  with  a  partially  adaptive  beamformer.  Two  of  the  articles  [20][19] 
discuss  a  method  based  on  the  use  of  the  GSC.  Both  of  these  methods  require  one  additional 
linear  constraint  for  both  narrowband  and  broadband  arrays.  Van  Veen  [43]  describes  a  method  to 
approximate  the  desired  response  which  uses  any  or  all  of  the  available  degrees  of  freedom. 

Partially  adaptive  beamformer  designs  are  analyzed  and  evaluated  in  several  sources.  Van 
Veen  [40]  mentions  two  design  methods  for  partially  adaptive  broadband  beamformers:  the  eigen- 
structure  based  and  the  power  minimization  based  design  procedures.  The  eigenstructure  based 
design  [42] [40]  problem  is  simpler  to  solve  than  the  power  minimization  problem,  but  requires  more 
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adaptive  weights  for  comparable  performance.  Power  minimization  based  design  [46][41][39]  gives 
the  best  interference  cancellation  performance  with  the  fewest  weights.  It  is  helpful  to  read  Liu  and 
Van  Veen’s  [26]  article  to  better  understand  the  power  minimization  based  design.  Van  Veen  [39] 
explains  that  beam  based  designs  are  not  justified  for  broadband  systems.  Partially  adaptive  beam¬ 
forming  is  out  cf  scope  in  this  thesis,  but  would  be  the  next  logical  extension  to  the  fully  adaptive 
case  examined  here. 

Van  Veen  [44]  wrote  about  soft  response  constraints  with  miuimum  variance  beamformers. 
LCMV  beamformers  use  hard  constraints  to  control  the  response  of  the  desired  signal.  If  the 
constraint  is  allowed  to  permit  distortion  of  the  desired  signal,  the  interference  cancellation  perfor¬ 
mance  can  be  improved.  Van  Veen’s  approach  is  to  minimize  the  variance  subject  to  a  constraint 
on  the  mean-square  response  (quadratic  constraint  on  the  weights). 

Several  articles  describe  coherent  and  correlated  narrowband  signals.  Li  and  Chang  [25] 
modeled  a  narrowband  multipath  signal  as  a  scaled  and  delayed  version  of  a  direct  path  signal.  The 
narrowband  multipath  scenario  generated  in  this  work  is  based  on  this  model.  Reddy  [34]  describes 
a  narrowband  model  for  two  coherent  sources.  Shan  [36]  and  Godara  [17]  also  cover  narrowband 
models  for  two  or  more  coherent  sources.  The  correlated  narrowband  signals  generated  in  this 
thesis  are  based  on  these  coherent  signal  models. 

Vural  compares  the  different  types  of  adaptive  array  processors  [47];  the  time  domain,  element 
space  processor  is  the  type  considered  in  this  study.  Adaptive  Technology,  Inc.  [29]  studied  the 
radome  multipath  reflection  problem.  Their  report  covers  the  problem  and  possible  solutions.  The 
idea  to  model  the  radome  as  two  flat  plates  originated  from  this  report. 
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III.  Preliminaries 


5.1  Introduction 

A  beamformer  weights  and  sums  data  from  an  array  of  sensor  elements  to  perform  spatial 
filtering.  Spatial  filtering  involves  separating  signals  which  have  similar  temporal  frequency  content 
but  arrive  from  different  directions.  Signals  are  classified  as  either  narrowband  or  broadband.  One 
definition  of  narrowband  is  in  terms  of  the  fractional  bandwidth.  The  fractional  bandwidth  is 
the  signal  bandwidth  as  a  percentage  of  the  carrier  frequency.  Signals  with  fractional  bandwidths 
much  less  than  1%  are  classified  as  narrowband  and  those  with  fractional  bandwidths  much  greater 
than  1%  are  broadband  [40].  Another  definition  of  narrowband  is  in  terms  of  the  time  bandwidth 
product  (TBWP).  The  TBWP  is  the  product  of  two  terms:  1)  the  temporal  aperture  or  the  time 
it  takes  a  wave  to  propagate  from  the  first  to  the  last  sensor  and  through  the  taps  and  2)  the 
bandwidth  of  the  signal.  Signals  with  a  TBWP  much  less  than  one  are  considered  narrowband 
and  those  with  a  TBWP  much  greater  than  one  are  broadband  [40].  A  narrowband  beamformer 
is  analogous  to  a  finite  impulse  response  (FIR)  filter.  While  the  response  of  a  FIR  filter  is  a 
function  of  frequency,  a  beamformer’s  response  is  a  function  of  both  angle  of  arrival  and  frequency. 
Narrowband  beamformers  are  much  easier  to  analyze  than  broadband  and  will  be  developed  first. 
The  broadband  beamformer  will  be  described  second.  The  beamformers  are  designed  according  to 
the  linearly  constrained,  minimum  variance  form. 

3.2  Basic  Definitions 

Vectors  and  matrices  are  used  throughout  this  thesis.  Vectors  are  denoted  by  bold  face  lower 
case  symbols,  and  are  assumed  to  be  column  vectors.  Matrices  are  denoted  by  upper  case  symbols. 
In  this  thesis,  A *,  AH  and  A _1  denote  the  operations  of  complex  conjugate,  Hermitian  or  complex 
conjugate  transpose,  and  inverse  of  a  matrix  A,  respectively. 
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3.3  Basic  Concepts 


3. 3.1  Narrowband  Beamforming.  Figure  3.1  ia  a  schematic  for  a  narrowband  beamformer. 
It  essentially  weights  and  sums  the  data  at  the  sensors  and,  accordingly,  is  sometimes  called  a 
weighted  delay-and-sum  beamformer  [11].  The  arrow  pointing  to  the  array  from  an  angle  0  with 
respect  to  the  array  broadside  (perpendicular  to  the  array)  represents  a  plane  wave  from  a  source 
in  the  far  field.  Plane  wave  propagation  implies  that  the  wave  amplitude  is  constant  at  all  points 


Figure  3.1  Narrowband  beamformer 


on  a  line  perpendicular  to  the  direction  of  travel.  The  beamformer  response  can  therefore  be 
characterized  only  by  delays  between  sensors,  and  not  amplitude.  The  elements  impart  a  sampling 
in  space  as  the  wave  propagates  through  the  array.  The  antenna  outputs  are  assumed  to  have 
been  sampled  in  time.  The  number  of  sensors  is  chosen  to  meet  the  Nyquist  criterion  for  spatial 
sampling,  which  implies 


d  < 


c 

2  •  /max 


•^min 
2  ’ 


(3.1) 


where  c  is  the  speed  of  propagation  in  the  medium  of  interest,  d  is  the  spacing  between  elements, 
/max  is  the  maximum  frequency  allowed  for  the  sources,  and  Amjn  is  the  wavelength  corresponding 


to  /max-  All  examples  in  this  thesis  use  normalized  frequency  as  expressed  by 


-»<w  <r  (3.2) 

Further,  c  and  d  are  often  normalized  to  one.  To  meet  the  Nyquist  criteria  for  space  and  time  we 
require  /max  =  1/2  or  Ami„  =  2. 

As  shown  in  Figure  3.1,  the  output  of  the  array  at  time  t  is 

N 

0(0  =  £«£*'{*),  (3.3) 

p=i 

where  N  is  the  number  of  antenna  elements.  A  compact  vector  form  of  this  equation  is 

0(0  =  w"x(0,  (3.4) 

where  the  data  at  the  sensors  are  represented  as 

*(0  =  M0  *a(0  "•  *Jv(01T>  (3-5) 


and  the  weight  vector  is 


w  =  [u>i  u>2  •••  u;w]T. 


(3.6) 


Assume  the  signal  prompt 'ng  toward  the  array  is  a  complex  sinusoid  with  direction  of  arrival 
(DOA)  0  and  radian  frequency  u.  The  zero  phase  reference  is  taken  to  be  the  center  of  the  array. 
This  point  need  not  coincide  with  a  physical  <*ntenna  element.  The  data  in  the  pth  antenna  channel 
due  to  a  unit  amplitude,  complex  sinusoid  in  the  far  field  is 


*p(0  =  e>"l‘+A,(*)l;  1  <p<N, 


(3.7) 


3-3 


where  Ap(0)  represents  the  additional  path  length  in  the  direction  of  propagation  from  the  aero 
phase  reference  to  the  pth  sensor: 


A ,(*)  =  ( 


N  +  l 
2 


~P) 


d  sin  5 
c 


(38) 


Substituting  this  into  equation  3.3  gives 


y(t)  = 

p= i 

N 

=  e** 

p=i 
j v 

p=i 

=  eiutr(6,w),  (3.9) 


where  r(0,  u)  is  the  beamformer  response.  The  response  can  be  expressed  in  vector  notation  as 


r(0,u>)  =  w*d(0,w). 


(3.10) 


Here  d(0 ,  w)  is  the  array  response  vector  defined  by 


d(0,W)  = 


1 

e -j'toi  J  (in  t 


(3-11) 
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By  defining  the  the  vector 


z"-' 

Equation  3.11  can  be  written  more  compactly  as 

d(0,w)  =  (i.iZ) 

From  equations  3.9  and  3.10,  the  beamformer  output  due  to  an  elemental  source  can  be  expressed 
as 

y(t)  =  e>“‘w*d(0,w).  (3.14) 

Here, 

x(O  =  e*-‘d(0,W),  (3.15) 

models  the  data  due  to  an  elemental  point  source  of  unit  amplitude. 

The  beampattern  is  defined  to  be  the  magnitude  squared  of  the  beamformer  response.  To  plot 
the  beampattern,  the  responses  at  the  angles  of  interest  are  needed.  The  row  vector  of  beamformer 
responses  can  be  written  as 

*■(">*)  =  [  r(u,e1)  r(U,e3)  ■■■  r(w,9N)  ] 

=  |  wBd(u,9i)  wHd(u/,0j)  •••  *rHd(w,0N)  j  •  (316) 
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Another  way  to  get  the  beamformer  response  is  to  take  the  discrete  Fourier  transform  (DFT) 


of  the  complex  conjugate  of  the  weight  vectors.  This  gives  responses  at  the  angles 


(317) 


To  see  this,  rewrite  equation  3.16  as 


r(w,  0)  = 


-  w# 


d(u,$i)  d(w,#2)  ••  d(w,0jv) 


(3.18) 


Substituting  the  array  response  vectors  into  equation  3.18  yields 


r(w,  0)  = 


_ 


1 


g-jwftin#,  in#j 


1 


c->ta»*«n*Ar 


e-;ui*(sr-l)«in«,  e- jui  i(N-l)  tin  l3  ...  e~i*i  (ff-l)tm  tN 


D,  (3.19) 


where  D  is  the  diagonal  matrix 


D  =  diag{e’"2(£[fi),in#I  e'ta’f(irA),in#a  ••• 


(3.20) 


The  matrix  in  the  center  of  equation  3.19  is  similar  to  the  DFT  matrix  which  is  given  by 


111  1 
1  e-iW#  •••  e-X1**"1)# 

1  e-W-D#  ... 


(3.21) 
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It  is  well  known  that  if  f  is  a  column  vector,  then 


gT  =  fTV  (3.22) 

is  the  DFT  of  fT.  Comparing  the  matrix  in  equation  3.19  with  the  matrix  representation  of  the 
DFT  in  equation  3.21,  we  see  that  when 

<9'=8in"1[^5^];  1<^N>  (3-23) 

then 

t(u,$)  =  whVD.  (3.24) 
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Figure  3.2  Broadband  beamformer 


S.S.S  Broadband  Beamforming.  Broadband  beamforming  is  similar  to  narrowband  except 
that  taps  are  added  to  allow  for  interference  cancellation  over  a  band  of  frequencies.  The  Nyquist 
criterion  for  spatial  sampling  is  met  as  in  the  narrowband  case.  The  delay  between  each  of  the  taps 
is  normalized  to  1  sec.  The  Nyquist  criterion  for  temporal  sampling  is 


2x 

=  —  >  2wmax- 


(3.25) 


With  T  normalized  to  one  second,  temporal  frequency  is  normalized  as 


—  *  <  u  <  x. 


(3.26) 
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Figure  3.2  is  the  schematic  for  a  broadband  beamformer.  With  the  broadband  beamformer,  there 
is  a  datum  at  each  of  taps.  The  equation  for  the  output  now  becomes 


K  L- 1 

y(0  =  E  E  wp,mxp(t'  m)>  (3-27) 

p=l  msO 

where  K  is  the  number  of  sensors  and  L  is  the  number  of  taps.  Delayed  versions  of  the  data  at 
each  sensor  can  be  written  as 

xp(t,  m)  =  eJ“f‘-mT+A»(*)],  (3.28) 


where  m  ranges  from  0  to  L  —  1  and  T  is  the  time  delay  between  adjacent  taps  illustrated  in 
Figure  3.2.  The  zero  phase  reference  is  at  the  center  of  the  array  as  in  the  narrowband  case. 
Substituting  the  above  equation  into  equation  3.27  gives 

y(t)  = 

p= 1  m=0 

=  ^  E  E  tn;,me-^[mT-A'(#)] 

p= 1  m=0 
K  l-l 

- 

p=  1  m=0 

=  e-jumT 

p=l  m=0 

=  e^'r^w).  (3.29) 


As  before,  r(0,w)  is  the  beamformer  response.  Again,  the  response  can  be  expressed  in  vector 
notation  as 

r(0,w)  =  w"d(0,u/),  (3.30) 
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where  d(0,  w)  is  now  expressed  as 


e-ju(l)iunt 


d(0,w)  =  eju'*<**l>,in# 


e-jtoi(X)Te-iui(l)*  sin  » 
c-iw(l)Te-;w(2)4»in# 

g-Jv(l)T g—jw(K)£  sin  • 


(3.31) 


c-iw(L-l)Tc-jw(l)*sin* 
e-ju/(L- l)Tc->ui(2)i  sin  « 

e-jw(L-l)Te-jw(K)f  tint 


and  the  weight  vector  is 


w  -  K.O  «>2,o  •••  u>k, o  toi(i  u^i  •••  ^,i  •••  wi,£-i  W2,l-\  •••  «>/r,£-i]T.  (3.32) 
This  array  response  vector  can  be  expressed  more  simply  as 

d(*,w)  =  *f*>-*,L(e-i-r)  0  (3.33) 
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where  the  symbol  ®  is  the  Kronecker  product  defined  in  [27].  An  example  of  a  Kronecker  product 


is 


A®  B  = 


ail B  ai2 B 
021 B  OL21B 


(3.34) 


where 


an  an 

A  = 

021  <»22 

and  B  is  a  matrix.  For  the  broadband  case,  we  evaluate  the  beamformer  response  as  a  function  of 
both  angle  and  frequency.  Therefore,  we  define  the  beamformer  response  matrix  as 


wHd(wi,0i)  wHd(wi,fl2) 
w"d(<j2,0i)  wffd(w2,02) 


w"d(wi,fljv) 

w"d(w2,0jv) 


w"d(wjy,0i)  mHd{uN,e2)  ■■■  w  Bd(uN,0N) 


(3.35) 


This  equation  may  be  solved  directly  for  the  response.  However,  as  in  the  narrowband  case,  an 
alternate  method  exists  for  calculating  the  response.  If  the  stacked  weight  vector  is  unstacked  and 
put  into  matrix  form,  the  DFT  of  the  rows  and  columns  can  be  used  to  obtain  the  beampattern. 
The  concept  is  the  same  as  in  the  narrowband  case  except  a  DFT  is  performed  on  each  row  (cor¬ 
responding  to  different  frequencies),  and  then  to  each  column  (corresponding  to  different  angles), 
of  the  beamformer  response  matrix.  This  case  is  more  complicated  than  narrowband  since  the 
angular  spacing  is  different  from  frequency  to  frequency.  At  each  frequency,  the  number  of  valid 
points  must  be  calculated. 
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3.4  Data  Generation 


All  the  data  used  in  this  thesis  are  simulated.  Since  the  data  are  assumed  WSS,  the  covariance 
matrix  of  x  can  be  written  as 

R,x  =  E[xx"].  (3.36) 

In  the  narrowband  case,  the  data  covariance  matrix  for  a  narrowband  signal  at  frequency  u>0  is 
represented  as 

Rxz  =  <T7d(0„Uo)dH(0„LJo)  +  Rnn,  (3.37) 

where  <r2  is  the  power  or  variance  of  the  signal,  0,  is  the  DOA  of  the  source,  w0  is  the  radian 

frequency  of  the  signal,  and  Run  is  the  covariance  matrix  of  the  noise.  When  the  signal  originates 

from  a  single  DOA  but  is  broadband,  the  expression  becomes 

Rzz  =  ~  /  s(w)d(0,w)d"(0,u-)  du  +  iz™,  (3.38) 

Jwen 

where  S(u>)  is  a  weighting  function  over  frequency,  called  the  power  spectral  density  (PSD),  and  Q 
is  the  frequency  extent  of  the  source.  The  most  general  expression  is 

Rxx  =  ^-  [  f  S(9,u)d{e.,w)dH{0„u)4wd0  +  Rnn,  (3.39) 

Jtte  Jwz n 

where  S(6,w)  is  a  weighting  over  frequency  and  DOA  and  is  also  a  PSD,  and  0  is  the  spatial 
extent  of  the  source.  The  noise  is  assumed  to  be  additive  and  white.  It  is  assumed  Gaussian  with 
covariance  <r2/,  where  <r2  is  the  white  noise  power  and  I  is  the  identity  matrix. 
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3.5  Linear  Constraint  Design 


The  simplest  type  of  constraint  is  that  of  a  point  constraint,  which  constrains  the  beamformer 
response  at  a  specified  frequency  and  DOA.  The  response  is  constrained  at  radian  frequency  w0 
and  angle  0,  as 

w"d(0,,u>,,)  =  r,  (3.40) 

where  r  is  the  desired  complex  response.  Several  linear  constraints  can  be  expressed  as 

C*  w  =  g,  (3.41) 

where  C  is  the  constraint  matrix  and  g  the  response  vector.  The  matrix  C  is  chosen  to  contain  M 
linear  constraints  so  the  above  equation  represents  M  linearly  independent  columns  in  N  unknowns. 
The  weight  vector  w  has  length  N,  where  in  the  broadband  case  N  =  K  x  L.  This  makes  the 
dimensions  of  C,  N  by  M .  If  N  =  M,  then  w  is  uniquely  determined  by  the  constraints.  To  ensure 
there  is  a  w  which  satisfies  the  constraints,  M  is  chosen  to  be  less  than  N.  An  example  of  multiple 
point  constraints  is 


[d(0»i.wo)  d(0,3,wo)  d(fl,s,w#)]ffw  =  [n  r2  r3]T.  (3.42) 


In  this  case,  M  =  3,  C  =  [d(0,i,wo)  d(0,2,wo)  d(0,3,w<,)],  and  g  =  [rx  r2  r3]T. 

A  more  efficient  type  of  constraint  in  the  broadband  case  is  the  eigenvector  constraint.  Eigen¬ 
vector  constraints  are  used  to  approximate  a  desired  broadband  response  using  only  a  few  linear 
constraints. 
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3.6  Linearly  Constrained  Minimum  Variance  (LCMV)  Beamforming 

3.6.1  LCMV  Criterion.  The  LCMV  criterion  applies  to  both  narrowband  and  broadband 
beamforming.  The  objective  is  to  pass  a  desired  signed  while  minimizing  power  from  undesired 
signals.  This  is  accomplished  by  constraining  the  beamformer  response  in  a  desired  direction  while 
minimizing  the  output  power  or  variance.  The  problem  to  be  solved  is 

min  P0  subject  to  CH w  =  g.  (3.43) 

W 

P0  is  the  expected  value  of  the  output  power  and  is  expressed  as 

PB  =  E[jy|2]  =  E[i/y*3  =  wHE[xx"]w.  (3.44) 

Thus,  the  problem  is  restated  as 

min  wffiZ*xw  subject  to  Ca w  =  g.  (3.45) 

To  begin  to  solve  this  problem,  guess  that  the  solution  is 

w,  =  R~lC[CH  (3.46) 

Equation  3.46  is  a  valid  guess  because  C  is  full  rank  and  Rax  is  positive  definite  (this  ensures 
the  inverses  exist)  [40].  Other  solutions  of  the  form  w  =  w,  +  a  meet  the  constraint  so  long  as 
CHn.  =  0.  Substituting  this  solution  back  in  the  objective  function  gives 

whRxxW  =  (wf  +  a")fl*x(w,  +  a) 

=  wf  fi**w,  +  wf  fl«a  +  a"  RxX  wf  +  aH  a 

=  wf  fl„w,  +  ga  \CH  R~xC\~lCH  Rgx  Rxx&  +  a  HRxxR^C[CaR^C]~1g  +  nHRxT* 
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-f  ««*w,  + 


=  w^H*rw#  +atfil„a,  (3-47) 

since  CHa  =  0.  We  see  that  all  other  possible  solutions  which  do  not  violate  the  constraint  make 
the  objective  function  larger.  Thus,  w,  gives  a  minimum. 

S. 6.2  Generalized  Sidelobe  Canceller.  The  generalized  sidelobe  canceller  (GSC)  is  another 
equivalent  form  of  the  LCMV  beamformer.  Figure  3.3  is  the  schematic  of  the  GSC  implementation 
of  the  LCMV  beamformer.  The  derivation  that  follows  is  from  Van  Veen  [40]. 


Figure  3.3  Generalized  sidelobe  canceller  schematic 

The  GSC  uses  an  orthogonal  basis  for  the  null  space  of  the  constraint  matrix  C.  Form  the 
matrix  C„  as  the  orthogonal  complement  to  the  space  spanned  by  the  columns  of  the  full  rank 
constraint  matrix  C.  Cn  is  size  N  by  N  -  M.  Together  the  columns  of  C  and  C„  span  the  entire 
N  dimensional  space.  The  weight  vector  can  be  expressed  in  terms  of  a  set  of  basis  vectors  as 

w  =  Cv  -  C„ w„,  (3.48) 
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where  v  is  Af  by  1  and  w„  is  N  —  M  by  1 .  Cv  and  C„w„  represent  the  components  of  w  in  the 
spaces  spanned  by  the  columns  of  C  and  Cn,  respectively.  Since  the  columns  of  C„  and  C  are 
orthogonal,  CHCn  =  0.  Applying  the  constraint  to  this  new  weight  vector  gives 

CH w  =  CHCy  -  CuCn w„  =  CHCv.  (3.49) 

Using  the  linear  constraint  equations  from  equation  3.43,  v  can  be  expressed  as 

v  =  (3.50) 

Substituting  this  back  in  equation  3.48  results  in 

w  =  C(CHC)-1g  —  C„w„.  (3.51) 

The  first  term  in  the  above  equation  only  depends  on  the  constraints  and  is  independent  of  the 
data.  This  part  is  implemented  in  the  upper  branch  of  the  GSC  and  represents  a  nonadaptive 
beamformer.  The  first  term  in  equation  3.51  is  sometimes  called  the  quiescent  weight  since  it  is 
the  part  independent  of  the  data.  It  is  represented  as 

w?  =  C(CHC)~1g.  (3.52) 

The  lower  branch  of  the  GSC  consists  of  an  unconstrained  adaptive  beamformer.  The  number 
of  adaptive  weights  is  .'V  ~  Af.  The  constraint  matrix  C  preserves  the  desired  signal.  Thus,  the 
matrix  C„  blocks  the  component  of  the  data  that  lies  in  the  space  spanned  by  the  columns  of  C. 
For  this  reason,  C„  is  referred  to  as  the  signal  blocking  matrix. 
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Using  the  GSC,  the  LCMV  problem  is  expressed  as  the  unconstrained  problem 


min  (wf  -  C„w„)H -  Cn wB). 


(3.53) 


The  solution  [40]  to  this  equation  is 


w„  =  {C*RmJ\)-lC*Rmm  w, . 


(3.54) 


The  reason  to  implement  this  form  is  that  the  direct  form  LCMV  beamformer  requires  a  constrained 
adaptive  algorithm,  which  is  more  complex  than  an  unconstrained  one. 
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IV.  Analysis  and  Simulations 


4-1  Introduction 

This  chapter  applies  concepts  from  chapter  3  to  the  simplified  RMR  problem.  Before  the 
actual  problem  is  addressed,  basic  narrowband  and  broadband  beamforming  concepts  are  tested 
and  compared  to  results  in  the  literature.  Since  the  implementation  of  the  guiding  equations  and 
algorithms  were  not  revealed  in  the  articles,  the  developments  are  included  here.  Much  time  will  be 
saved  by  using  these  derivations  and  algorithms  in  the  future.  The  actual  algorithms  are  written 
using  Matlab.  The  files  are  included  in  the  appendix. 

4-2  Narrowband  Bcamforming 
4-2.1  Data  Generation. 

4-2. 1.1  Uncorrelated  Sources.  Complex  data  are  generated  for  the  narrowband  case. 
The  covariance  of  a  single  narrowband  source  is  represented  [24]  as 

R.  =  <r*  d(0„u.)d{0„u,)a.  (4.1) 

Here  the  data  x  are  generated  as  a  random  amplitude  times  the  array  response  vector.  For  two 
uncorrelated  sources  this  becomes 


x  =  s  +  n 


—  ajdi  +  ajdj  +  n, 


(4.2) 
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where  s  is  the  received  signal  vector,  the  o/s  are  random  signal  amplitudes,  and  the  dj’s  are  array 
response  vectors.  This  equation  can  be  represented  as 


x  =  Bm  +  n,  (4.3) 

where  a  =  [ai  a2]T,  B  =  [di  dj],  and  n  is  uncorrelated  noise  represented  as  [i»i  n2]r.  The  vector 
a  is  sometimes  called  the  signal-in-space  vector  [24].  Assume  ai  and  a2  are  statistically  independent 
and  the  vectors  a  and  n  are  independent  and  normally  distributed  with  zero  mean.  Further,  n  is 
Gaussian  with  covariance  <r2I.  The  data  covariance  matrix  can  be  computed  as 

=  E[xx*] 

=  E[(Ba  +  n)(aHBH  +  n")] 

=  E[Baa*Bw]  +  Efna^B*]  +  E[BanH]  +  E[nn"] 

=  BE[aaH]BH  +  B„» 

—  BRtaB^  +  Ann 

=  A„  +  Rnn  i  (^4) 

where  R„  is  the  spatial  correlation  matrix  of  the  signal  and  Rnn  =  <t„In  with  In  the  identity 
matrix  of  size  N.  The  matrix  Raa  is  the  spatial  correlation  matrix  of  the  signal-in-space  vector  and 
is  expressed  as 


R 


aa 


E[«a*] 


E 


a2aj 
a2a*  a2a2 


E[aiaJ]  EfaiaJ] 
E[a2oJ]  E[a2Oj] 
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of  0 
0  of 

Pi  0 
0  Pi 


where  E[a<a?]  =  E[|«* I3]  =  of  =  Pi  Notice  that  R,a  can  be  expressed  as 


R.M 


0  of 
[di  d,] 

l' 

df 

[«fdi  ofd ,] 

<*? 

of  didf  +  ofd2df 


of  0 

d? 

0  of 

df 

+  ftin 


+  ftiii 


=  f2,i  +  R,i  +  Rnn  ■ 


(4.5) 


(4.6) 


This  is  an  expected  result  for  uncorrelated  sources. 

4. 2. 1.2  Correlated  Sources.  Godara  [17]  represents  two  correlated  signals  as 


s  =  >/Pi0idi  +  VP3(Pe°i  +  \A  —  |P«l*a3)d2,  (4.7) 

where  <n  and  02  are  independent  and  have  zero  mean  and  unit  variance.  Also,  pi  and  pj  are  the 
powers  of  signals  one  and  two,  respectively.  In  the  correlated  model,  the  off  diagonal  elements  of 
Raa  are  proportional  to  the  correlation  coefficient.  The  correlation  coefficient  [37]  is  defined  as 


Pe  V^IaiaflEIa^l 


(4.8) 
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This  can  be  rearranged  as 


ElajoJ]  =  pcy/pipi  .  (4.9) 

The  matrix  Raa  is  Hermitian  symmetric  so 

E[a2aJ]  =  ply/piPi  ■  (4.10) 


Now,  compute  the  data  covariance  matrix  with  the  Rta  as  above: 


Rii  —  [di  dj] 


Pi  Pcy/PlPi 
Ply/pipi  Pa 


d? 

d? 

4"  ftin 


=  [pidi  +  p*«Wpipj  pedi>/PiP3  +  <*aPa] 


d? 

d? 


4"  Rnn 


=  Pididf  4-  P^dj^/PiPadf  4-  Pcd^pipjdf  4-  pjdjd?  4-  Rnn- 


(4.11) 


4-2- 1.8  Multipath  Model.  The  multipath  model  used  simulates  a  specular,  single 
bounce  reflection  off  a  flat  plate  into  the  antenna  array.  Figure  4.1  shows  the  set  up.  The  plate 
on  the  right  does  not  reflect  energy;  signals  pass  through  unaffected.  Signals  from  positive  DOAs 
reflect  off  the  left  plate  according  to 


PrV)  =  (|}a«(*).  (412) 

where  pr(9)  is  the  reflection  coefficient,  u(0)  is  the  unit  step  function,  and  $  is  the  angle  from  the 
broad  side  of  the  array  in  radians.  This  correlation  coefficient  function  of  6  was  chosen  because  it 
approximates  that  of  an  actual  radome.  When  signals  arrive  at  broad  side  (perpendicular  to  the 
array),  there  are  few  or  no  reflections.  Signals  reflect  more  as  the  angle  from  broad  side  increases 
(100  %  reflection  at  90°  from  broad  side).  The  reflection  process  is  assumed  to  only  attenuate  the 
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nroph—  mtw ncsunsyc— ■ 

Figure  4.1  Two  plate  multipath  model 

signal  after  it  bounces  off  the  plate.  No  random  phase  is  introduced  by  the  reflection  surface;  it  is 
assumed  to  be  perfectly  smooth.  The  wave  starts  out  planar  and  remains  that  way  after  it  bounces 
off  the  wall.  The  reflected  or  indirect  path  wave  is  modeled  as  a  scaled  and  delayed  version  of  the 
direct  path  wave.  For  a  narrowband  signal,  the  signal  part  of  the  data  observed  by  the  sensors  is 
given  by 


8  =  Oidi  +  02^2 

=  aidi  +  praie-;*d2 

=  aidi  +  PrOie-' ^Afld2,  (4.13) 

where  A  R  is  the  path  length  difference  between  the  direct  path  and  indirect  path  wave  at  the  zero 
phase  reference  of  the  array. 
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The  signal  covariance  matrix  for  this  model  is 


R„  =  E[ss*] 

=  E[(ojdi  +  we-**- ^da)(a;df  +  pW*A*d? )] 

=  EtoiaJJdxdf  +  Elaiajy^^pf  drdf  +  E^Ke-' *A*dadf  +  E[a1a;]|pr|2dJd? 
=  ffj'didf  +  <rJe'^fAi*p;didf +  (Tjpre"''^Ai,d3dj/ +  <ri|pr|ad2df.  (4.14) 


The  data  covariance  matrix  for  the  above  signal  covariance  matrix  can  be  written  in  vector  notation 


as 


Rxx  =  [di  dj] 

This  gives  Raa  as 


a\  e\p'rti*** 

df 

alpre-jtf~AR  <r\\pr\2 

d? 

+  Rnn- 


Raa  — 


a\  <r\p*^R 

<r?pre~HfA*  <r%\p,\2 


(4.15) 


(4.16) 
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Figure  4.2  Narrowband  beamformer  with  eight  interferere;  SNR  =  INR  =  60  dB 

4-2.2  Point  Constraints.  Point  constraints  are  used  with  the  narrowband  beamformers. 
In  all  cases  the  beamformer  response  is  constrained  to  unity  in  the  direction  of  the  desired  signal. 
Figure  4.2  displays  the  beampattern  for  a  beamformer  constrained  at  the  DOA  arcsin(O.l)  with  a 
signal  environment  consisting  of  eight  uncorrelated  interferere  from  direction  sines  of  -0.85,  -0.55, 
-0.25,  0.05,  0.25,  0.45,  0.65,  and  0.85  in  uncorrelated  noise.  The  frequency  of  all  signals  is  1/2  (or 
u  =  x).  The  interferere  have  power  levels  60  db  relative  to  the  noise.  Figure  4.2  shows  this  result, 
which  is  the  same  as  that  obtained  by  Van  Veen  [40].  The  file  bfnb9s.m  contains  the  code  for  this 
problem.  Using  the  GSC  form  of  the  LCMV  beamformer  provides  the  same  results  as  the  direct 
form. 

The  GSC  form  with  point  constraints  is  used  for  the  multipath  simulation.  The  frequency 
of  all  signals  is  again  1/2  (w  =  x).  The  multipath  case  is  set  up  with  the  interferer  arriving 
from  +45  degrees  and  reflecting  into  the  desired  look  direction,  which  is  chosen  as  -45  degrees. 
This  is  the  worst  case  for  the  multipath  in  that  all  the  reflected  signal  will  reflect  into  the  main 
beam.  When  the  interference  signal  reflects  into  the  main  beam,  a  null  formed  in  the  DOA  of  the 
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Figure  4.3  Narrowband  beampattern  for  multipath  case 

interferer  will  not  solve  the  problem.  The  beamformer  must  cancel  the  direct  path  signal  with  the 
indirect  path  signal.  The  more  correlated  these  two  signals  are,  the  better  they  can  be  cancelled. 
Since  the  multipath  reflection  signal  is  a  scaled  and  delayed  replica  of  the  direct  path  signal,  the 
two  are  correlated.  The  indirect  interference  signal  can  be  destructively  added  to  the  direct  path 
interference  signal  because  they  are  correlated.  Since  the  reflection  coefficient  at  45  degrees  is  1/4 
(from  equation  4.12),  the  reflected  signal  is  one  fourth  the  amplitude  of  the  direct  signal  and  delayed 
in  time.  To  cancel  these  signals,  the  beamformer  forms  a  beam  four  times  as  large  in  amplitude 
in  the  direction  of  the  reflected  signal  as  in  the  direct  path  direction  (see  figure  4.3).  Since  this 
is  a  beampattern  plot  (magnitude  squared),  the  difference  in  the  beam  gains  is  12  dB.  To  see  the 
amount  of  cancellation,  the  output  power  of  the  beamformer  is  calculated.  This  gives  the  response 
of  the  entire  system  (including  plates)  to  unit  amplitude  far  field  sources.  Figure  4.4  is  the  output 
power  of  the  beamformer  with  the  quiescent  weights,  while  figure  4.5  is  the  output  power  with 
the  adapted  weights.  The  Matlab  files  outputpownbq.m  and  outputpownbw.m  contain  the  code  to 
compute  and  plot  the  output  power.  Figure  4.5  shows  that  the  interference  signal  at  45  degrees 
has  been  cancelled. 
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DOA 


Figure  4.4  Output  power  with  quiescent  weights;  SNR  =  60  dB 


Figure  4.5  Output  power  with  adapted  weights 
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4-3  Broadband  Beamforming 


4-3.1  Data  Generation. 

4-3. 1.1  Uncorrelated  Sources.  Real  data  are  generated  for  the  broadband  ease. 
Constant  PSDs  are  assumed  for  all  sources;  this  weights  all  frequencies  evenly  and  is  the  worst  case 
assumption  [40].  The  PSD  used  is 


S(u)  =  .](")  +  ^X[w.,u,k](w), 


where  x(w)  is  the  indicator  function  defined  as 


X[*,,](«) 


1  if  x  <u  <y, 

< 

0  otherwise. 


As  defined  earlier,  the  array  response  vector  is 


(4.17) 


(4.18) 


d(w,0)  =  e'-^* ™,*L(e-iwT)  ®  *K(e->ut “*). 


(4.19) 


Signal  covariance  matrices  for  the  uncorrelated  sources  with  no  spatial  extent  are  calculated 
as  follows: 

R$  =  ~  /  S(w)d(w,^)d/f(w,5)dw  +  ~  f  S(u>)d(u/,6)da  (u,0)dw 

=  7"  /  U'  ^,6)dH(u,e)du>+^-  H  d(u,0)dH(u,0)du 

1  f  —  U/  * 

4?r  J 

+  j-  /  [^^J»in,*i(e-><*'T)®*K(e-','ta'^in,)][c;'u'£fi'“n#*t(e">‘,,3’)®»K(e_ito’*‘in#)]Hdu; 

Jut9 
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(4.20) 


The  two  terms  in  the  Kronecker  product  can  be  expressed  as 


*L(e->u'T)9%(e-iuT)  = 


1  eid^T  ...  gjCL-i^T 

e-j(iVT  j  ...  giCi-aVT 

e-j(t-i)«T  ...  2 


=  A(w)  (4.21) 


and 


*K(e~M"ntf)*2(e-i"',in#)  _ 


c-i(l)u-*«nS 

g—j{K— 1)«  *  tin  t 

B(u,9). 


ej(l)"f  «n* 
1 


The  signal  covariance  matrix  can  now  be  compactly  written  as 


e>(Jf-2)w*»i  in# 


1 


(4.22) 


1  i  /■»» 

R>  =  -^J  a(u)  ®  B{u,  9)  du>+—J  A(w)  ®  S(«,  0)  du».  (4.23) 


The  above  matrix  can  be  expressed  as 


R,= 


4x 


f-ut  °oB(u,  9)  dw  /_  "*  Oi B(u,  9)  dw 
XT"*  a-iB(u,  9)  dw  /j;;  aoB(u,  9)  du> 

«-(i-i)*(«,*)  dw 


/_Jk‘  a(L-i)S(w,0)  dw 

r“*a(£-2)S(u;,tf)da; 


JZm  a0B(u>,9)  6u, 
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£  a0B(u,,6)  du  C.^B(u,0)du  •••  £  a(t.-l)B(w,$)  dw 

1  £';a_1B(W,fl)dW  £  aoB{u,0)  du  •••  j£  a(L.7)B(w,0)  dw 


(4.24) 


fZ*  a-(L-i)B(u,0)  dw 


r>B(u,«)  d« 


where  a<,  =  e^“wT.  The  above  matrix  is  block  Toeplitz.  The  entire  matrix  can  be  formed  by 
calculating  the  blocks  of  the  first  row  and  column  and  creating  a  Toeplitz  matrix.  This  is  done 
using  Matlab.  nr  <e  block  on  the  i[h  principle  diagonal  of  the  above  matrix  is 


R(ii)  =  /  "*  ^uTB(u,0)  dw  +  r  e>i'wTB(u,0)  dw.  (4.25) 

J  —  Ulh 


The  elements  on  the  tj*  principle  diagonal  of  R(ii)  are 


r(i1(i2)=  /  “‘cii*"V‘1"«"n,dw+^-  F'  e’i'"Teii**i"n,dw. 

J -oil 


(4.26) 


The  variable  ij  ranges  from  -(L-l)  to  (L- 1)  in  increments  of  1.  The  variable  i2  ranges  from  -(K- 1) 
to  (/f-1)  in  increments  of  1.  Each  element  of  the  signal  covariance  matrix  is  in  closed  form  as 


r(n,i2)  =  -L  [  ^  e3AnT+i3i.in«}du  +  ±_  r‘^[..T+<1*.in«ldw 

4ir  J_W(1  4x  yw. 

e-i«.l<i-T+»a  *.in*l  _  e-ju>*[»i  T+»a  *»in»l  _j_  g;w»[»i  T+ia*«in#]  _  gjw.hi  T+ta  Jiin*] 

4xj[«i  •  T  +  »2  •  |  sin  6} 

8in[w6r(»i  T+i 2  •  |  sin  5)]  -  sintu/oT^ii  •  T  +  «2  •  f  sin  5)] 

2x(ii  •  T  +  »2  •  |  sin  0) 

_  wjT  •  sinc[w6(»i  •  T  +  »2  •  |  sin  5)]  —  wor  •  sinc[w„(»i  •  T  +  *2  •  §  sin  0)] 


(4.27) 


with  the  sine  function  defined  as 


sinc(z)  = 


sin(z) 


(4.28) 


4-12 


Equation  4.27  is  used  to  compute  necessary  elements  of  the  signal  covariance  matrix.  The  file 
mscmnew.m  implements  this  equation. 

Van  Veen  [39]  described  another  method  of  computing  the  signal  covariance  matrix  for  two 
uncorrelated  signals  at  a  single  location.  First,  he  defined  R,  for  a  single  source  as 

R,  =  /  pJ(w)d(w)dH(w)do;)  (4.29) 

Jn 

where  p2(w)  is  the  source  PSD  and  fi  is  the  source  frequency  extent.  He  approximates  the  above 
integral  by  a  Riemann  sum  and  avoids  integration  in  the  following  equation: 

j 

R,  »  ^  p2  (w)d(w)dH  (w)  Am 
i= i 

w  AT3Ah,  (4.30) 

where 

A  =  [d(wi)  d(w2)  •  •  •  d(wj)] 

T2  =  diag{p2(w1),p2(u;2),--,p2(w/)}Aw.  (4.31) 

Now,  the  two  uncorrelated  sources  and  noise  can  be  expressed  as 

T?  0  Ah 

Rtt  « [Ai  A2)  1  +Rnn.  (4.32) 

o  rl  a f 
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This  equation  is  used  for  negative  and  positive  frequencies  to  obtain  real  data.  Note  that  d(-w)  = 
d*(w)  and  d"(-w)  =  dff*(w)  =  dT(w).  Reed  data  are  obtained  with  the  following  equation: 


Rxx  «  [^i  ^2] 


Ea 

2 


A? 

A? 


+  [a\a-2  } 


? 


r3 

La 

2 


AT 

AT 


+  Rnn- 


(4.33) 


The  file  RJntur.m  computes  the  signed  covariance  matrix  using  this  method. 

Next,  suppose  the  source  has  spatial  and  frequency  extent.  The  PSD  function  used  is 


S{u>9)  = 


(4.34) 


In  this  case,  the  covariance  matrix  involves  a  double  integral  and  is  defined  by 


‘  S(w,  0)d(w,  0)dH(u,  0)dud0  +  jf  *  S(w,  0)d(u,  0)dH(u,  0)du>d0 


1  /•*  n 


=  4* I  »J  * e)iude d(w’ *)d"  (w- 

=  ^  J  J  W*[eJ“^^“,*£®«K][^w^'"n'*i®«K]"du;d<> 

4T  Je.  Jw, 

=  ^  J  *'  J  “'[*L*?]®[*K*5)dvdO+~  [9  K9HK]dud0,  (4.35) 


where  the  arguments  of  9l  and  9  k  have  been  dropped  and  are  e~’wT  and  e~>wc respectively. 
This  can  be  solved  as  in  section  4.3. 1.1.  Each  element  of  the  signal  covariance  matrix  is  in  closed 
form  as 


r(n ,  i2)  =  L  J  '  J  W*  e>“[<lT+<a i  *Jdu/d0  +  jf *  c*-P*T+*»« “  *5du>d0 
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=/ 


e-jw.[>i  T+»J  * ain #]  _  --jw»[i,  T+ia  Jtin#]  [»*  Juji,  T+i,  *m»]  _  .;««[«,  T+»a  f  tin  #] 

- - d0+  /  - - d0 

-it  4 T  +  ij  -4  gin  fl]  7#.  4x;[i1  T  +  t2  •  7  «n  tf] 


-z 

-j:. 


-4*j'[t'i  •  T  -  ij  •  4  sin  0] 

*»  e-jw.[«i  T-ia  * tin#]  _  e~jwt[ii  T-ta  gtin#) 


*•  e-jw.[i,  T-»a  *  tin#]  _  *->«»(»,  T-«a  f  tin  #]  r*i  -jwalii  T+«a  *  tin  #]  _  .jw.lii  T+«a  *tm#] 

dff  +  /  . . . . . ■—■■  dfl 

7#.  4*il*i  T +  ij  fsintf] 


- s - dff  +  (  - ^ - ——j - ^ - d  e 

4jtj[*i  •  T  -  ij  •  f  sin  0]  7#.  4*j[ti  •  T  +  »j  •  f  sin  0] 


*»  gjwjii  T+«a  f  tin#)  _  giu.tii  T+ia  *tin#] 


r* t  T-ij  f  tin*]  _  e-jw»J«i  T-»a  *  tin»J  gjvtlii  T+17  £  tin*]  _  gjw.lti  T+ta  j  tin#] 

/  l  rn  •  d  ill  *  —  -f-  'T'  ,  -  d  /ll 

7# . 


4*j[«i  ■  T  -  *2  •  7  sin  0] 


4*7 [ii  •  T  +  »2  •  7  sin  0] 


7^  e->oi«[ij  T-»a  *tin#J  _  e-ju»k(ii  T -«a*  tin #]  gjuiahi  T+«3  *  tin #]  _  gjw.tii  T+i,  f  tin#] 

4*j[ii  •  T  -  ia  ■  7  sin  0]  4*j[»j  •  T  + 12  •  7  sin  0]  ^ 


i=i 


(4.36) 

The  above  equation  is  implemented  in  Rdb.m.  The  file  Rtdb.m  calls  the  function  Rdb.m  and 
computes  the  entire  covariance  matrix. 
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Figure  4.6  Two  plate  multipath  model 


4.8. 1.2  Multipath  Model.  Figure  4.6  shows  the  geometry  for  the  broadband  multi- 
path  model.  It  is  the  same  set-up  as  in  section  4.2. 1.3  for  the  narrowband  case.  The  path  length 
difference  between  the  direct  path  ray  and  the  indirect  path  ray  at  each  array  element  is  calculated 
to  determine  the  array  response  to  the  reflected  signal.  The  DOA  0  is  the  angle  from  the  broad 
side  of  the  array  as  depicted  in  figure  4.6.  The  path  length  difference  is  computed  based  on  the 
geometry  of  figure  4.6: 

Pi  +  P3  =  2(P  +  i  •  d)  sin  0,  (4.37) 

where  i  is  the  array  element  number  with  respect  to  the  phase  center  of  the  array  as  defined  in 
figure  4.6,  and  d  is  the  array  sensor  spacing.  The  path  lengths  P\  and  Pj  sum  to  give  the  total 
path  length  difference  between  the  direct  and  indirect  path  rays  at  the  ith  sensor.  With  this  path 
length  difference  at  each  element,  the  reflected  signal  can  be  represented  as 


where  c  is  the  propagation  speed.  The  total  signal  can  be  represented  as  the  sum  of  two  signals, 
the  direct  path  signal  and  the  reflected  one  as  follows: 


S  =  S1  +  s3. 


(4.39) 


The  data  covariance  matrix  can  now  be  computed  as 


R%  —  E[s*H]  +  Rnn 

—  E((Sl  +  Sj)(s®  +  •?))  +  Rnn 

—  E^isf]  +  E[siS]f]  +  E[«jsf]  +  E[s2s^]  +  R„n 

—  fill  +  fi»!ia  +  fi»j»  i  +  fii,  +  Rnn  i  (4.40) 

where 

R%,  =  f  S(w)di(w,0)df(w,0)dw 

Jn 

ft.,.,  =  /  S(w)di(w,0)<lf(w,0)dw 

Jn 

ft.,.,  =  /  S(w)d3(w,  0)d?  (u,8)Au 

Jn 

R„  =  f  S(w)dj(w10)djr(u/,0)dw 

Jn 

Rnn  =  all,  (4.41) 


with  d\  the  array  response  vector  of  the  direct  path  signal  as  originally  defined,  d3  the  array 
response  vector  of  the  reflected  signal  modified  according  to  equation  4.33,  and  <r*  the  noise  power. 
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The  array  response  vectors  are  expressed  as 


d,  =  *(e->“T)® 


1 

g-/lw  J  tin  t 


and 


d2  =  *(e"->“T)®pr 


g/ltoi*  (in  tin  # 

1  .e-i*  af  »u>* 

e->  lw*  (in  #e-jw  sin  # 


=  '9(c~*wT)  ®  Pr 


gjlu*  till  nn  #e-ju>2*  sin  * 

1  .  e->w2f  «n# 

e->  Ini  J  «in  #£-jh>2£  (in  tfjvl  £  (in  # 


_  #(c-j<-3’)®pre->'-2?"n# 


e->lw*«in# 

1 

gjlwftin# 


(4.42) 
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The  last  equation  can  be  written  as 


da  =  *.e->“3?  •**<!!(«,-*). 


(4.43) 


The  above  equations  for  the  array  response  vectors  indicate  that  the  direct  path  signal  covariance 
matrix  is  the  same  as  previously  defined.  The  covariance  matrix  of  the  second  (reflected)  signal 
is  just  p*  times  the  previously  defined  signal  covariance  matrix  computed  for  the  negative  of  the 
angle  of  the  first  signal.  The  expressions  for  the  cross  term  signal  covariance  matrices  are  new  and 
must  be  determined.  Substituting  the  PSD  and  second  array  response  vector  in  4.41  leads  to 


R. 


1*3 


=  ^Pr  j  ^  “%(«,  0)d“(w,  -0)du  +  ~Pr  f  *dj(u/,  0)d?(w,  -0)du> 

^  J 

+ ±-Pr  f*  e>"*« ' ““»*L  ®  9)()Hdu 

=  7~Pr  f  W‘  e’u2S™t[e>u*rli™,*L®*K][eiwE*li,ia'*L  ®  ***]<**' 

J 

+Tpr  ®  **"]dw 

=  ^PrJ_ 

+^Pr  ^ •t?*K'l*™t[9L9f)  ®  [9k*Jc) dw, 


(4.44) 


where  9l9%  is  defined  in  equation  4.21.  Again,  the  arguments  of  9i  and  9k  have  been  dropped 
and  are  e~>wT  and  c-i"f  “n,J  respectively.  Define  D(u,0)  as 


D(w,0)  =  ^K-1)i^t9K9TK  = 


^(/f-lluifiinS  ej(A'-2)a>**in#  ...  j 

gj(K-2)uinn«  gj(/r-3)u>*«  n#  ...  g-jlwfun# 

2  e-;'ltoi*»in#  ...  e-j(JC-l)u>*  •“>* 


(4.45) 
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The  cross  covariance  matrix  can  then  be  written  compactly  as 


=  ±-prJ  eiui$™l  A(w)®D(u>,$)dw+±pr  JU'  e’u2$™*  A{u)<2>D(u,,$)<L>.  (4.46) 


Ai,!,  can  now  be  written  as 


4t 


Pr 


e>“3?*in#a0£>  dw  eJ“3? “%£>  dw 

e^w3?,",#a_1i?  dw  /"“*  e>“2^"n#a0D  dw 


e>“3f  •in*a.(L_1)Z>  dw 


^^*a{L.7)D  <l> 

/_■“*  e>"*$™ta0D  du> 


+  4^ 


j£  giwJf  .in  #ao£>  ^  jr*  e^3*  “  %£)  da> 


14 


4  *in  • 


rv 


a(i_i)Z)  dw 


fa  e*u7j‘"BtaoD  dw  •  J£  e>“2?  dw 


£  Aw 


J^4e»"3f""#a0Dda; 


(4.47) 


where  =  e-,,,wT  and  the  arguments  of  D  have  been  dropped.  The  above  matrix  is  block  Toeplitz. 
The  whole  matrix  can  be  formed  by  calculating  the  blocks  of  the  first  row  and  column  and  creating 
a  Toeplitz  matrix.  This  is  done  using  Matlab.  The  block  on  the  i[h  principal  diagonal  of  the  above 
matrix  is  a  Hankel  matrix  of  the  form 

R(h)  =  /  *  <J“3?“n  V,'*"T£>(w,  A)  Aw  +  P  <J‘*'3^'in  V,'1“’TI>(w,f>)  dw.  (4.48) 
The  elements  on  the  minor  diagonal  of  R(ii)  are 


r(»i ,  h)  =  f  '  ju2$™ite’ixwTe>i*,*'in9AM+  f  ,ei»Jf««»eji,wTtii,"J“*du.  (4.49) 

•'“Wb  v  W« 
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The  variable  t'i  ranges  from  -(L-l)  to  (L-l)  in  increments  of  1.  The  variable  12  ranges  from  -(AM) 
to  (/f-1)  in  increments  of  1.  Each  element  of  the  signal  covariance  matrix  is  of  the  following  form: 


— pr  /  *  eju>{iiT+2f«n#+ia!«m#]<iw  +  J_  ef«[«,T+2f  .in  #+<a  *  *m 

4*  4t 

e-jw.[<i'T+2£  sin  f+ia- *«■>*]  _  g-jw.ft,  T+2£  tmt+iji  (in  #] 


T  •  T  -t-  sin  0  +  ij  ■  |  sin  0] 

gjwtpi  T+2£.in»+»a  *.in»]  _  gjw.li,  T+2^  un  #+«j  *  sin  #] 

P'  4xj[ii  •  T  +  2*  sin  0  + t'2  •*  sin 0] 

sin[wiT(»i  •  T  +  2^  sin  0  +  1*3  •  &  sin  0)] 

Pr  2ir(»i  •  T  +  2^sin0  -f  tj  •  |  sin0) 

sin[w0T(»i  •  T  +  2^  sin  0  +  12  •  ^  sin  0)] 

Pr  2x(ii  •  T  +  2^  sin  0  +  12  •  |  sin  0) 

u/jT  •  sinc[w»(*i  •  T  +  2^  sin  0  +  12  •  7  sin  0)1 
p  - ; - s - £ - ” 

Pr  2  T 

w0T  ■  sinc[wa(ii  •  T  +  2^  sin  0  +  ij  ■  £  sin  0)] 

-Pr  2t  ■ 


(4.50) 


The  file  mscml2h.m  computes  this  cross  signal  covariance  matrix.  The  other  cross  covariance 
matrix,  R,2„ ,  is  calculated  in  the  same  manner.  The  structure  of  is  the  same  as  Its 

elements  can  be  expressed  as 


, .  .  .  UbT  •  8inc[wj(«i  •  T  -  2^  sin  0  +  »2  •  7  sin 0)1 
r(*i ,  t2)  =  pr - 27 - - 


-  pr- 


uaT  ■  sinc[wa(ii  T  -  2y  sin  0  +  »2  •  §  sin  0)] 


2* 


The  file  mscm2Ih.m  computes  this  cross  signal  covariance  matrix. 


(4.51) 
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4-3.2  Eigenvector  Constraints. 


4-S.2.1  General  Case.  The  broadband  constraint 

wH  d(w,0)  =  r(u,0)  for  <  w  <  1*4,  6a  <6  <6^  (4-52) 


consists  of  an  infinite  number  of  constraints  as  opposed  to  a  point  constraint  in  the  narrowband 
case.  Writing  many  constraints  at  specified  frequencies  and  angles  would  be  inefficient  because  each 
constraint  uses  up  one  degree  of  freedom.  Also,  the  constraints  would  hold  at  only  the  specified 
frequency  and  angle  pairs.  There  would  be  no  control  over  the  response  at  frequency  and  angle 
pairs  other  than  the  specified  ones.  For  this  reason,  Van  Veen  [40]  defines  eigenvector  constraints 
which  minimize  the  total  squared  error  between  the  desired  response  and  the  actual  response  over 
the  entire  frequency  and  angle  extent  while  using  up  a  fixed  number  of  degrees  of  freedom.  Defining 
rd  as  the  desired  response,  the  total  squared  error  between  the  desired  and  actual  response  is 


-hCC  | rs(0,u)  —  w/fd(0,u'j'l2dud0  +  J  j"  | rj(0,u)  —  wKd(0,w)|2dwdfl. 

(4.53) 


The  problem  is  to  minimize  this  total  squared  error  subject  to  keeping  a  fixed  number  of  constraints. 
The  factor  £  is  included  for  convenience.  With  the  previous  definition  of  R, . 


R,  =  /  *  /  *  d(0,w)dH(0,w)dud0  +  -1  ['*  r*  d(0,u)dH(0,w)dud0, 

J-»t  J-ub  JUm 


the  problem  is  stated  as 


min  [wTi?,w  -  —  wT  /  /  rj(0,w)d(0,w)du>d0 

w  4ir  J_tk  J_wt 


1_ 

4x 


ri(0 ,  w)d*T  (0,  w)dwd0 
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-  s*Trr  r*(0,w)d,T(0tw)dwd0^  w],  (4.54) 

where  the  transpose  operator  is  substituted  for  the  Hermitian  since  the  data  and  weights  are  real 
and  the  terms  excluding  w  are  eliminated.  After  taking  the  gradient  with  respect  to  w  and  setting 
the  equation  equal  to  zero,  this  can  be  arranged  as 


1  fib  fU/b 

R,w  =  —  J  j  rj(-9,-w)d(-O,-L>)dbid0 


+  J_  /'*  r‘r<J(-e,-w)d*(-»,-«)dwdO 
8*  J» .  Ju. 

+  —  f  r  rj(8,w)d(0,u>)dwdO  +  f  f  i\i(0,w)d*(0,w)dwd0. 
8*  Ji.  Ju.  Jt.  Ju,. 


(4.55) 


Define  f  as  follows: 


i  rib  ruib  i  rtb  ruib 

f=F-/  /  »*3(-0,-u/)d(-0,-u/)dwd  0+i-  /  r<j(— 0,  — w)d*(— 0,  — w)dwd0 

8""  Jt.  Ju,.  8tt  Jft  JWt 


tb  ruib 


1  rib  ru, b  |  rVb  fU>b 

+  —  /  /  r3(0,w)d(0,w)dwd0  +  —  /  /  rd(0,  w)d*  (0,  w)dwd0.  (4. 

8*  Jea  Ju,  8x  Jem  Ju. 

By  properly  choosing  the  desired  response,  equation  4.56  can  be  solved.  Constrain  the  response  to 
have  unit  gain  and  linear  phase  with  the  response  defined  as 


Ob  fWb 


56) 


t&  =  e-X^T 


(4.57) 


This  amount  of  delay  is  required  for  a  natural  beamformer  response.  If  you  do  not  allow  enough 
delay  in  the  constrained  response,  the  beamformer  will  not  have  time  for  weighting  and  summing 
enough  of  the  delayed  data  at  the  sensors  and  taps  to  meet  the  constraint. 
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Each  of  the  four  integrals  is  solved  separately  to  obtain  f .  The  first  integral  is  solved  by 


f  r  r*d{-9,  -w)d(-0,  -w)dw 
Jt€  Juim 


1 

^01(^5^)*  iin# 

■  //‘-',VVT 

«in# 

-IJ 


L 

J 

L 

c-j(±riVr 

eM^)*  «n* 

c2UI('£yi_K+l)f  .in* 

dwd  0 


6uA6. 


Each  of  the  elements  in  the  above  vector  can  now  be  expressed  as 


7(ii)  =  J  J  e~iilwT 

ejw(£f±)*  tint 

.in* 

dwd  6  = 

W> 

W-o 

eJw(^tjJ._JC+l)*.in* 

7(«ja_K+i) 

where  each  of  the  7_,as  is  of  the  form 

7-< a  =  J  J  e-iilwTe~*i3W''in*<iwd0 


(4.58) 


(4.59) 


(4.60) 


with  »i  fixed.  The  variable  «i  ranges  from  to  —  L  +  1.  The  variable  »2  ranges  from  —K^1 
to  ~~K£~  +  K  —  1.  The  other  integrals  can  be  solved  in  a  similar  manner  to  the  first  one.  The 
resulting  equations  and  ranges  of  indices  are  provided  below.  The  elements  of  the  second  integral 
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in  equation  4.56  are  of  the  form 


7,-a  =  j  e>i'wTe’i*‘i™,du 


(4.61) 


with  «i  fixed.  The  variable  «'i  ranges  from  to  —  L  +  1.  The  variable  »2  ranges  from 
to  +  K  -  1.  The  elements  of  the  third  integral  in  equation  4.56  are  of  the  form 


7<a  =  J  e’iluTe><aUilitttdv 


(4.62) 


with  ii  fixed.  The  variable  ij  ranges  from  to  ^y  —  L  +  1.  The  variable  t2  ranges  from 
to  —  K  +  1.  The  elements  of  the  fourth  integral  in  equation  4.56  are  of  the  form 


7-ia 


(4.63) 


with  *i  fixed.  The  variable  «i  ranges  from  ^y  to  -  L  +  1.  The  variable  «2  ranges  from  ^y 
tofiyi-AT+l. 

Since  I  had  already  written  the  following  algorithm,  the  first  plus  second  integrals  were 
computed  and  the  third  plus  the  fourth  to  find  the  elements  of  f : 


/(»i ,  i2)  =  j*k  er“[*'.T+.a4.in»]dw^  d0  +  _L  j'  e-MhT+i,i .in'ld^  dO 


1  /■**  e>"»[nl'+*3f  «n#]  _  er<*,«[*i3’+«a**in»]  j  e-jw*[«iT+ia£  nn«]  _  e-j««[»iT+ia * tin*] 


=  -/ 

Sir  J». 


l  r* 

j[tiT  +  t2f  sinS]  8x  J,m  -j[*iT  +  *2£sin0] 

1  /**  ei<*’t[>iT+»a J»in I]  _  gjw.tnT+iji.in#]  c-;w»[»iT+ia * «in S]  _  e-2w.[«iT+«3f  tin#] 


d0 


=  —  f 


j[iiT  +  »2£sin0] 


-i[»iT  +  i2^sin0] 


d0 


1  »in  »)  _  e7w.[«iT+i3  Jiin  »]  e-r«»(»,T+<a*  tin  »]  _  e-ju>«[i,T+«a*  tin  #] 

«  -  )  - - - 4.  - - -  AS. 

8x  j[ixT  +  i2 1  sin  0]  -j[iiT  +  i2f  sin  S] 


(4.64) 
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The  above  equation  can  be  used  to  compute  the  first  pair  and  the  second  pair  of  integrals  since  the 
pairs  have  common  ranges  for  the  indices  of  ii  and  «2.  Each  of  the  pairs  of  integrals  is  of  the  form 


f (M ,  **)  =  r-  r  r‘e>to'[<,T+<l«“n#1dwdfl+  /  ‘  r\->“ li*T+<’*“n#Jdwd0.  (4.65) 

Jt.  Ju.  ®T  Jt.  Ju>. 


The  files  ftdb.m  and  fdb.m  implement  the  above  integrals  with  a  Riemann  sum.  Equation  4.55  can 
be  written  in  vector  notation  as 

R.  w  =  f.  (4.66) 

This  is  in  the  same  form  as  the  linear  constraints  previously  defined  as 


Ca  w  =  g, 


(4.67) 


with  R,  in  the  place  of  the  Hermitian  of  the  constraint  matrix  and  f  in  place  of  g.  Perfor  <  ing  a 
low  rank  approximation  on  R,  will  be  more  efficient  by  using  less  degrees  of  freedom  than  all  of 
the  columns  of  R,.  Keeping  most  of  the  energy  of  the  eigenvalues  provides  a  good  approximation 
to  higher  rank  matrix,  R,.  The  signal  covariance  matrix  can  be  written 

R,  *inirVH,  (4.68) 

where  Er  is  a  diagonal  matrix  containing  all  or  most  of  the  nonzero  eigenvalues.  Based  on  Buckley’s 
[5]  time  bandwidth  product  for  sources  at  a  single  angle,  a  formula  is  developed  which  provides 
good  results  for  the  multiple  angle  case.  The  number  of  constraints,  D,  used  in  the  low  rank 
approximation  is  computed  as 

D  =  \[(K  -  1){  -  sin  e2  -  -  sin  *  }  +  (L  —  l)T|2(/2  +  (4.69) 

c  c 
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where  fx]  is  defined  as  the  integer  greater  than  or  equal  to  x.  Substituting  equation  4.68  in  equation 
4.66  and  noting  that  U  is  a  unitary  matrix  (UHU  =  /)  yields 


UZrVH  w 
UhU'LtVh* 

VHw 

VHw 


=  f 
=  UBt 
= 

=  f'. 


(4.70) 


This  in  again  in  the  form  of  the  constraint  equation 


CH  w  =  g, 


(4.71) 


with  C  =  V  and  g  =  f '  =  L~lUHf.  These  values  can  be  substituted  in  the  optimum  weight 
equation  to  determine  the  steady  state  weights. 

When  the  angle  of  the  source  is  a  fixed  point,  f  is  defined  as 


f  =  ^  f'  r'd{60,  -u,d(0Ol  -u)dw  rd(0o,  -w)d*(0o,  -u)Au, 


^  f'  r2(««,,w)<i(S0)w)du;  +  p  rd(0o ,  w)d* (tf0 ,  w)dw ,  (4.72) 


in  equation  4.66.  In  this  case,  the  number  of  constraints  used,  D,  is  based  on  capturing  99.99 
percent  or  more  of  the  energy  represented  by  the  eigenvalues  of  R,.  The  fraction  of  energy  in  the 
eigenvalues  is  calculated  as 


fx  = 


(4.73) 


4-27 


where  N  =  K  x  L.  The  elements  of  f  are  found  as 

/(*i,«a)  =  T^-{w*sinc[w»(tiT  +  ij-  sin0)]  -  uasinc[w.(>iT  +  *j-  sintf)]}.  (4  74) 
in  c  c 

Equation  4.74  holds  for  the  multipath  model  for  sources  having  frequency  extent  with  a  slight 
modification.  The  desired  response  must  be  modified  to  account  for  the  additional  delays  through 
the  beamformer  associated  with  the  reflected  (delayed)  signal  into  the  array.  The  beamformer  needs 
time  to  weight  and  sum  the  signals  which  are  delayed  relative  to  the  direct  path  signal,  which  is 
the  only  signal  accounted  for  in  the  previous  definition  of  rj.  All  of  the  delayed  signals  need  to  be 
combined  to  form  the  output  of  the  beamformer.  If  all  of  the  delayed  data  at  the  sensors  and  taps 
are  not  combined,  the  constraint  response  cannot  be  met.  The  modified  response  is 

r<l  =  e--»(47A)“Te--»‘*’«*in#.  (4.75) 

Computing  the  same  four  integrals  above  with  the  modified  response  yields 

1  P  d 

f(«i ,  ij)  =  —  {w»sinc[wi(ii7’  +  —  •  sin  9  -4-  ij  -  sin  0)] 

P  d 

—  w08inc[wa(«iT'  H —  •  sin0  +  ij-sinfl)]}.  (4.76) 

c  c 

The  file  myfcenref.m  implements  the  eigenvector  constraints  for  the  multipath  case. 


4-28 


4-3.2. 2  Correlated  Sources /Muhtpath  Environment.  The  correlated  and  multipath 
signal  environments  are  described  next.  The  signals  of  the  first  case  are  perfectly  correlated  (pc  =  1). 
This  example  was  performed  to  evaluate  the  cancellation  performance  for  two  correlated  undesired 
signals  when  one  of  them  arrives  in  the  direction  of  the  main  beam.  The  two  undesired  signals 
are  assumed  to  be  uncorrelated  with  the  desired  signal.  The  file  bfbbrefn.m  contains  the  code  for 
this  problem.  The  array  has  16  elements  equally  spaced  along  a  line  at  one-half  the  wavelength 
of  the  highest  frequency  employed,  x.  The  tap  spacing  is  normalized  to  one  second.  This  again 
normalizes  frequency  on  the  interval  — zr  to  x  radians.  Twelve  eigenvector  constraints  are  employed 
to  ensure  unit  gain  and  linear  phase  over  the  frequency  interval  [2x/5, 4x/5]  at  negative  45  degrees 
(main  beam).  The  undesired  signals  are  located  at  positive  45  degrees  (right  of  array  broad  side) 
and  negative  45  degrees  (left  of  broad  side)  with  power  levels  of  20  db  relative  to  the  white  noise 
power  level.  Figure  4.7  shows  that  a  beam  is  formed  in  the  direction  of  the  signal  at  +45  degrees 


-1  -05  0  03  1 

•moi  oca 

Figure  4.7  Gain  of  beamformer  with  correlated  signals  as  a  function  of  DOA  for  constraint  over 
frequency  interval  [2x/5, 4x/5]  at  -45  degrees  plotted  at  nine  frequencies 

with  the  same  amplitude  level  as  the  signal  at  -45  degrees.  The  beamformer’s  weights  adapt  to 
minimize  the  output  power.  When  the  two  undesired  signals  are  perfectly  correlated,  a  null  is  not 
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formed.  The  signals  are  cancelled  by  destructive  addition.  In  this  example,  the  responses  to  the 
two  undesired  signals  are  formed  so  they  will  cancel  (approximately  the  same  amplitude  but  close 
to  160  degrees  out  of  phase).  The  cancellation  is  computed  using  the  following  equation  for  the 
beamformer  output  power: 

P0  =  ww/2riw.  (4.77) 

The  output  power  for  this  case  is  -43.19  dB,  which  shows  that  the  signals  have  cancelled  well. 
This  cancellation  eliminates  most  of  the  undesired  signal  energy  at  the  beamformer  output.  In 
the  case  of  a  radar  in  the  presence  of  two  correlated  jammers  (each  uncorrelated  with  the  desired 
signal),  this  cancellation  effect  degrades  the  effectiveness  of  the  jamming.  If  the  signals  happen  to 
be  correlated  desired  sign&ls  (as  can  occur  with  near  field  multipath  scattering  into  the  array),  the 
desired  signal  will  cancel.  This  is  highly  undesirable. 

The  last  example  is  for  the  case  of  the  radome  simulation.  A  signal  arrives  from  +45  degrees 
and  bounces  off  a  flat  plat  and  into  the  array  at  -45  degrees.  The  direct  path  signal  is  assumed 
to  be  unaffected  by  the  radome.  The  metal  plate  attenuates  the  signal  but  is  assumed  to  impart 
no  random  phase  to  the  signal.  The  reflected  signal  is  a  scaled  and  delayed  replica  of  the  direct 
path  signal.  This  simulates  a  signal  bouncing  into  the  main  beam.  The  file  bfbbrefnn.m  contains 
the  code  for  this  problem.  The  array  has  16  elements  with  25  taps  per  element.  The  SNR  is  20 
dB.  The  reflected  signal  is  attenuated  by  one  fourth  after  reflecting  from  the  plate.  Ten  eigen¬ 
vector  constraints  are  employed  to  ensure  unit  gain  and  linear  phase  over  the  frequency  interval 
[3*/5, 4x/5]  in  the  desired  signal  direction  of  -45  degrees.  Figure  4.8  shows  the  unity  constraint  at 
-45  degrees  and  the  beam  formed  at  +45  degrees  at  one  fourth  the  power.  This  plot  shows  a  twelve 
dB  drop  because  it  is  a  plot  of  the  beamformer  response  squared.  The  interference  cancellation  is 
demonstrated  in  figures  4.9  and  4.10.  Figure  4.9  is  the  output  power  of  the  beamformer  with  the 
quiescent  weights  while  figure  4.10  is  the  output  power  with  the  adapted  weights  due  to  the  source 
arriving  from  +45  degrees  in  addition  to  the  reflection  at  -45  degrees.  These  plots  are  generated 


Figure  4.10  Normalized  output  power  with  adapted  weights 

in  the  same  way  as  the  beampattems.  In  either  case,  the  weights  are  first  computed  for  the  envi¬ 
ronment  of  interest.  For  the  beampattern,  the  responses  due  to  a  narrowband  source  swept  across 
the  directions  of  arrival  are  plotted  versus  DOA.  Each  of  these  beampattern  plots  is  for  a  specific 
frequency.  In  the  power  plots,  the  output  powers  due  to  a  broadband  source  swept  across  the 
directions  of  arrival  are  plotted.  They  include  the  effects  of  the  simulated  radome.  For  this  reason, 
the  power  plots  are  like  wide  band  responses.  The  Matlab  files  outputpowq.m  and  outputpoww.m 
contain  the  code  to  compute  and  plot  the  output  power.  Figure  4.9  can  be  though  of  as  installed 
antenna  pattern  of  a  nonadaptive  beamformer  where  the  large  sidelobes  at  the  negative  of  the  main 
beam  angle  are  due  to  RMR.  Figure  4.10  can  then  be  thought  of  as  an  effective  installed  antenna 
pattern  of  an  adaptive  beamformer.  The  actual  installed  beampattern  with  adapted  weights  (fig¬ 
ure  4.8)  shows  that  the  beamformer  actually  forms  a  beam  in  the  interferer  direct  path  direction 
of  +45  degrees.  Figure  4.11  is  a  plot  of  the  beamformer  output  power  as  a  function  of  the  ratio  of 
the  temporal  duration  of  the  source  to  the  delay  between  the  direct  and  indirect  path  at  the  phase 
center.  Cancellation  is  affected  by  both  the  temporal  duration  and  the  delay.  If  the  delay  between 
the  different  paths  is  small  enough  the  correlation  can  approach  unity.  The  delay  was  previously 
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Figure  4.11  Cancellation  as  a  function  of  the  ratio  of  temporal  aperture  (T(0))  to  path  delay 
(del(0)) 

defined  as 

del{0)  =  2Psin(0).  (4.78) 

The  delay  gets  smaller  with  smaller  angles  (closer  to  broadside).  Cancellation  improves  as  the 
correlation  increases.  For  example,  the  output  power  with  L  =  2  taps  is  -16.93  dB  at  0  =  45°  and 
-21.15  dB  at  0  =  15°.  For  a  given  angle,  the  temporal  duration  increases  with  increasing  number 
of  taps  (for  a  fixed  number  of  sensors).  Cancellation  performance  improves  with  more  taps.  For 
example,  at  6  =  15°,  the  output  power  is  -21.15  dB  with  2  taps  and  -35.43  dB  with  15  taps. 

Another  parameter  which  affects  the  cancellation  performance  is  the  source  frequency  band¬ 
width.  With  the  difference  in  path  lengths  so  small,  the  major  contributor  to  decorrelation  of  the 
signals  is  the  wide  bandwidth  of  the  interference  signal.  The  wave  field  coherence  decreases  with 
bandwidth  [22].  This  is  due  to  the  lack  of  coherence  at  time  differences  of 

At  >  I  (4.79) 
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where  B  is  the  frequency  bandwidth.  For  example,  at  6  =  75°  with  two  taps,  the  output  power  is 
-16.7  dB  for  a  source  frequency  bandwidth  of  0.1  and  -9.43  dB  for  a  source  bandwidth  of  0.2. 
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V.  Conclusions  and  Recommendations 


The  results  obtained  show  that  if  the  interferer  signal  reflecting  off  the  radome  into  the 
main  beam  is  well  correlated  with  the  direct  path,  the  signal  can  be  cancelled  without  forming  a 
null  in  the  direction  of  arrival.  With  the  narrowband  multipath  environment  (zero  bandwidth), 
the  output  power  plot  shows  the  cancellation  is  about  50  dB.  With  the  broadband  multipath 
environment  (normalized  frequency  interval  of  [3ir/5,4ir/5]  radians),  the  output  power  plot  shows 
the  cancellation  is  about  20  dB.  The  cancellation  performance  improves  significantly  as  the  number 
of  taps  is  increased.  For  the  multipath  environment  generated  in  the  previous  chapter  with  a 
SNR  =  20  dB,  with  DOAs  of  ±20°  instead  of  ±45°,  the  output  power  is  -21.15  dB  with  2  taps 
and  -35.43  dB  with  15  taps.  As  the  signal  frequency  extent  increases,  the  degree  of  cancellation 
decreases.  For  the  same  multipath  environment  as  noted  above  with  DOAs  of  ±75°  and  two  taps, 
the  output  power  is  -16.7  dB  for  a  source  frequency  bandwidth  of  0.1  and  -9.43  dB  for  a  source 
bandwidth  of  0.2.  Depending  on  the  frequency  bandwidth,  a  moderate  number  of  taps  can  be 
added  to  obtain  good  cancellation  performance. 

Desired  signal  energy  will  cancel  due  to  RMR  sidelobes  in  a  manner  similar  to  interference 
signals.  In  this  study,  concentration  was  placed  on  cancelling  reflected  interference  signals  while 
constraining  the  beamformer  response  in  the  direction  of  desired  signals.  The  data  were  generated 
to  consist  of  interferers.  In  an  actual  radome,  desired  signals  may  reflect  into  RMR  sidelobes.  Since 
these  reflected  signals  may  be  correlated  with  the  one  in  the  desired  direction  (main  lobe),  part 
of  the  desired  signal  might  cancel.  The  amount  of  cancellation  will  be  dependent  on  the  actual 
reflection  coefficient  of  the  radome,  the  bandwidth,  and  the  directions  of  arrivals  of  the  signals. 

The  results  obtained  so  far  for  the  GSC  are  encouraging.  The  beamformers  employed  in 
this  thesis  are  fully  adaptive.  These  beamformers  utilize  all  of  the  available  adaptive  degrees  of 
freedom,  and  are  difficult  to  implement.  Each  adaptive  degree  of  freedom  typically  requires  a 
receiver  channel  These  channels,  in  addition  to  the  channel  needed  for  the  fixed  beamformer,  are 
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not  easily  added  to  an  airborne  radar.  The  next  logical  step  to  this  thesis  is  to  evaluate  interference 
cancellation  performance  using  only  a  small  number  of  the  adaptive  degrees  of  freedom.  This  is  the 
subject  of  partially  adaptive  beamforming  and  is  discussed  in  many  of  the  reference  articles.  Van 
Veen  [40]  has  written  several  articles  on  partially  adaptive  beamforming. 

The  method  of  partially  adaptive  LCMV  beamforming  does  not  require  pre-selection  of  the 
particular  elements  which  comprise  the  auxiliary  antennas  as  in  the  case  of  the  multiple  sidelobe 
canceller.  The  MSC  assumes  a  known  main  and  auxiliary  channel  designation.  In  a  GSC,  a 
transformation  is  performed  on  the  N  —  M  available  adaptive  degrees  of  freedom  to  reduce  the 
number  of  adaptive  weights.  A  likely  interference  scenario  is  developed  which  contains  information 
such  as  the  maximum  number  of  interferers,  their  locations  and  frequency  ranges,  white  noise  level, 
etc..  This  information  is  used  to  design  the  transformation  matrix,  which  is  built  into  the  array. 
Since  the  locations  of  the  jammers  are  not  known,  the  locations  can  be  guessed  to  range  over  the 
angular  field  of  view  of  the  array  (-90°to  +  90°).  The  beamformer  design  based  on  these  locations 
should  be  evaluated  to  see  how  well  it  works  for  different  test  scenarios. 

Another  area  which  should  be  investigated  is  beampattern  control.  The  normal  beampatterns 
obtained  from  the  LCMV  method  have  typical  -13  dB  down  sidelobes.  Most  radar  applications 
require  much  lower  sidelobe  levels.  If  a  partially  adaptive  beamformer  is  designed,  the  unused 
adaptive  degrees  of  freedom  can  be  used  to  approximate  a  desired  quiescent  response  [20]  [43].  The 
quiescent  pattern  can  then  be  approximated  to  beampatterns  typically  obtained  by  deterministic, 
nonadaptive  arrays. 

In  this  thesis,  all  frequency  is  in  terms  of  normalized  frequency.  Intermediate  frequencies 
(IFs),  the  IF  bandwidth,  and  the  IF  sampling  rate  of  a  typical  X-band  radar  should  be  used  next. 

The  radome  models  can  be  expanded  to  include  more  complex  elliptical  models  which  more 
closely  resemble  the  actual  radomes.  These  models  which  do  not  make  use  of  plane  waves  are  very 
difficult  to  develop.  Actual  array  data  at  the  individual  elements,  if  available,  should  be  used. 
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Appendix  A.  MATLAB  Code 


All  Matlab  code  used  to  support  this  thesis  was  written  by  the  author. 

Xbf .response. m 
clear 

a  s  20;  %  number  of  sensors 
lambda  =  2;  X  wavlength 
d  =  1;  X  interelement  spacing 

v  =  (l/n)*«xp(-j*2*pi*d/lambda*.l).“[0:n-l). ';  X  weight  vector  at 
X  sin (DO A)  =0.1 

theta  =  pi* [-1/2: .00006:1/2] ;  X  range  of  angles  to  compute  response  at 
gamma  =  ezp(-j*2*pi*d/lambda*sin(theta)) ;  X  form  of  array  response 
X  vector  element 

D  =  (ones(n, l)*gamma) .~([0:n-l] . '*ones(l,max(size(gamma))));  X  array 
X  response  vectors 
X 

X  add  a  Hamming  taper  to  the  beamformer  weights 
Xk  =  -10:9; 

Xwh  »  .64  +  .46*cos(pi*k/10) ; 

Xw_tap  =  wh . ’ . *w ; 

Xr  =  w_tao’*D; 

X 

r  =  w’*D;  X  responses 
plot(sin(theta) . 10*logl0(abs(r) .“2)) 
grid 

Xplotnb.m 

X  n  is  the  number  of  sensors,  I  is  the  number  of  angles  minus  one, 

X  d  is  the  interelement  spacing,  lam  is  the  wavelength,  and  w  is  the 
X  weight  vector 

function  [R]  =  plotnb(n,l,d,lam,w) 
theta  =  pi/I* (0:1) -pi/2; 

D  =  (ones(n,l)*(exp(-j*2*pi*d/lam*sin(theta))))."([0:n-l]  . ’*ones( 1,1+1)) ; 
R  =  w'*D; 

plot(sin(theta) , 10*logl0(abs(R) . ‘2)) 
grid 

zlabeK'sin  of  D0A’) 
y label ('Gain  in  dB’) 

Xplotfftnb.m 

X  I  is  the  number  of  points  in  the  FFT  and  w  is  the  weight  vector 
function  [R]  =  plotfftnb(I.w) 

R  =  fftshift(fft(conj(w) ,■)) ;  X  I  point  FFT  of  conjugate  of  weight  vector 
R(I+1)  =  R(l);  X  add  one  point  to  response  to  enable  plotting  from  -1 
X  to  +1 

k  =  C-l:l/(I/2) :1] ;  %  range  of  values  for  sin(theta) 
plot(k,10*logl0(abs(R) .*2))  %  plot  gain  of  response  vs  sin(theta) 
grid 

zlabeK’sin  of  D0A’) 
y label ( 'Gain  in  dB') 
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Xbinb9t.il 

clear 

n  *  20;  X  number  of  eluents 
Xb  =  200;  X  number  oi  data  vectors 
laabda  =  2;  X  wavelength 

Xsigaa  =  le-3;  X  standard  deviation  oi  the  noise 
sigBa_sq  =  le-6;  X  variance  or  power  oi  noise 

R_aa  =  eye(9) ;  X  spatial  correlation  matrix  oi  signal-in-space  vector 
X 

X  decomposition  oi  R_aa 
XCv.e]  =  eig(R.aa); 

XA  =  v*sqrt(e)*v'; 

X 

XI  =  sigaa*randn(n,B) ;  X  m  noise  vectors 
d  *  1;  X  intereleaent  spacing 

X 

X  array  response  vectors  oi  signal  and  interierers 
ds  =  exp(-j*2*pi*d/lambda*(0. 1))  ([0:n-i].’); 

dil  =  exp(-j*2*pi*d/lanbda*(-.85))  .*  ([0:n-l].’); 
di2  =  exp(-j*2*pi*d/lambda*(-.65))  ([0:n-l] . ’) ; 

di3  =  exp(-j*2*pi*d/laBbda*(-.25))  .*  ([0:n-l].’); 
di4  =  exp(- j *2*pi*d/lambda*(0 . 05) )  ([0:n-l].’); 

di5  =  exp(- j*2*pi*d/lambda*(0 . 26) )  ([0:n-l].'); 

di6  =  exp(-j*2*pi*d/lambda*(0.45))  .*  ([0:n-l].'); 
di7  *  exp(-j*2*pi*d/laabda*(0.66))  ([0:n-i].’); 

di8  *  exp(-j*2*pi*d/lambda*(0.85))  .*  ([0:n-l].'); 

X 

c  =  ds;  X  constraint 
i  =  1;  X  desired  response 

w_o  =  c*inv(c’*c)*f ;  X  quiesent  weight  vector 
X 

X  make  C_n  orthogonal  to  c  and  w_o 
cc  =  [c  w_o] ; 

[O.S.V]  =  svd(cc) ; 

C_n  =  U(: ,2:n); 

C_n’*c 

C_n’*w_o 

X 

H  =  [ds  dil  di2  di3  di4  diS  di6  di7  di8] ;  X  direction  matrix 

XX  -  H*A*randn(9,m)  +  I;  X  m  data  vectors 

XR_xx_hat  =  l/m*X*X' ;  X  sample  covariance  matrix 

R_xx  =  H*R_aa*H’  +  sigma_sq*eye(n) ;  X  data  covariance  matrix 

Xnorm(R_xx_hat  -  R_xx) 

RI  =  inv(R_xx);  X  inverse  oi  R.xx 

w  =  RI*c*inv(c’*RI*c)*i;  X  optimum  weight  vector 

w_n  =  inv(C_n’*R_xx*C_n)*C_n’*R_xx*w_o;  X  adaptive  weights 

w2  =  w_o  -  C_n*w_n;  X  weight  vector  in  GSC  iorm 

XRIh  =  inv(R_xx_hat) ; 

Xwh  =  RIh*c*inv(c’*RIh*c)*f ; 
iignre(l) 
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Xplot f f tnb2 ( 1024 ,w,wh) ; 
plot?  t !tnb(  1024 , *) ; 

?ignro(2) 

plotnb(n, 10000, d, lambda, w) ; 
figureO) 

plotnb(n, 10000, d, lambda, *2); 

?igur«(l) 

Xplotfftwb.m 

X  K  is  the  number  of  sensors,  L  is  the  number  of  taps  per  sensor,  I  is 
X  the  number  of  points  in  the  2-dimensional  FFT,  f  is  the  frequency, 

X  and  u  is  the  weight  vector 
function  [R]  =  plotf f tab (K , L , I , ? , w) 

V  s  reshape (a, K, L) ;  X  unstack  the  weight  vector  in  a  matrix  where  the 
X  columns  are  delayed  array  output  vectors 

R  =  fftshift(fft2(conj(W),I,D);  X2-dinensional  FFT  of  conjugate  of 
X  weight  matrix 
X 

X  add  another  row  and  column  to  allow  plotting  sin(theta)  from  -1  to 
X  +1  and  normalized  radian  frequency  from  -pi  to  +pi  radians 
R(: ,1+1)  =  R(: ,1); 

R(I+1 , : )  =  R(l, :); 

X 

X  determine  how  many  points  to  use  based  on  the  frequency 

fn  3  round(I*(f+.5)+l);  X  frequency  number 

Is  »  (I+l)-(fn+l)+l; 

incs  =  floor((fn-(ls+l)+l)/2); 

k  =  [-l:l/incs:l] ; 

X 

plot(k,10*logl0(abs(R(ls+l:fn,fn)) .*2))  X  poll  out  column 
X  corresponding  to  frequency  number  and  rows  corresponding  to  good 
X  spatial  points 
grid 

xlabelC 'sin  of  D0A') 
ylabel(’Gain  in  dB’) 

Xmscmnew.m 

X  K  is  the  number  of  sensors,  L  is  the  number  of  taps  per  sensor,  the 
X  is  the  angle  the  source  is  arriving  from,  fl  and  f2  are  the  lower 
X  and  upper  limits  of  the  normalized  frequency  of  the  source,  doc  is 
X  the  interelement  spacing  divided  by  the  speed  of  propagation 
function  [R]  =  mscmnev(K,L, the, fl, f2.doc) 
wb  3  2*pi*f2;  X  normalized  radian  frequency 
wa  =  2*pi*fl;  X  normalized  radian  frequency 
th  3  the*pi/180;  X  angle  in  radians 
x  =  doc*sin(th); 

R  3  zeros (K*L,K*L) ;  X  fora  data  covariance  matrix  with  all  zeros 
for  il  =  -L+1:L-1;  X  for  each  value  of  il,  compute  row  and  column 
for  i2  =  0:K-1;  X  find  row  and  column  of  matrix  for  each  il 
row(l,i2+l)  =  (l/(2*pi))*(wb*msinc(wb*(il+i2*x))-wa*msinc(wa*(il+i2*x))) ; 
col(l,i2+l)  =  (l/(2*pi))*(wb*msinc(wb*(il-i2*x))-wa*msinc(wa*(il-i2*x))) ; 
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«nd 

Rsub  «  toeplitz(col .row) ;  Xform  toeplitz  matrix 
if  il  =*  0; 

R  *  R  +  kron(eye(L) , Rsub) ;  X  add  matrices  to  main  diagonal 
•lsaif  il  >  0; 

X  add  matrices  to  diagonal  above  main  diagonal 
R  =  R  +  kron(diag(ones(l,L-il)lil),Rsnb); 
elseif  il  <  0; 

X  add  matrices  to  diagonal  below  main  diagonal 
R  =  R  +  kron(diag(ones(l,L+il),il),Rsub); 
end 
end 

XR.intur.m 

X  li  is  the  number  of  frequency  increments,  thl  is  angle  of  first 
X  source,  th2  is  the  angle  of  the  second  source,  fl  and  f2  are  the 
X  lover  and  upper  frequencies  of  the  sources,  sigl.sq  and  sig2_sq  are 
X  the  variances  or  posers  of  the  first  and  second  sources 
function  [RI]  =  R_intur(K,L,Ii,thl,th2,fl,f2,sigl_sq,sig2_sq) 
omegal  =  2*pi*fl; 
omega2  =  2*pi*f 2; 

inc  =  (omega2  -  omegal)/li;  X  find  frequency  increment 
omega  -  omegal  +  inc*(0:li);  X  frequency  points 
thetal  =  thl+pi/180; 
theta2  =  th2+pi/180; 

dl  «  zeros(K*L,Ii+l);  X  set  array  response  vectors  of  source  1  to  zeros 
d2  =  zeros (K+L, Ii+1 ) ;  X  set  array  response  vectors  of  source  1  to  zeros 
X 

X  for  each  frequency  compute  the  array  response  vector  of  the  sources 
for  il  =  l:Ii+l; 

dl(: ,il)=kron(exp(-j*omega(il)).*([0:L-l]  ■ ’),  ... 
exp(j*omega(il)*sin(thetal)) ."([(K-l)/2:-l: ((K-l)/2)-K+l]  . ’)) ; 
d2( : ,il)=kron(exp(-j*omega(il)).“([0:L-l] . •  •  • 
exp(j*omega(il)*sin(theta2)).*([(K-l)/2:-l:  ((K-l)/2)-K+l] .  ’)); 
end 
X 

A1  =  dl; 

A2  =  d2; 

Gam_sq  a  0.6*eye(Ii+l)*inc;  X  flat  PSD  for  real  data 
R.aa  a  [Gam_sq  zeros(Ii+l, Ii+1) ;zeros(Ii+l, Ii+1)  Gam.sq] ; 

X  covariance  matrix  over  positive  frequency  for  two  uncorrelated  sources 
Rip  =  [sigl_sq*Al  sig2_sq*A2] *R_aa* [A1 ' ;  A2’]; 

X  covariance  matrix  over  negative  frequency  for  two  uncorrelated  sources 
Rim  =  [sigl_sq*conj(Al)  sig2_sq*conj (A2)]*R_aa*[Al. • ;  A2.*3; 

RI  =  Rip  +  Rim;  X  covariance  matrix  for  two  uncorrelated  sources 


XRdb.a 

X  oa  and  ob  are  the  lover  and  upper  radian  frequencies,  thl  and  th2 
X  are  the  lover  and  upper  angles  of  the  source,  li  is  the  number  of 
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X  angle  increments,  il  and  i2  are  indicies  from  HTdb.n 
function  [R]  =  Rdb(oa,ob,thl,th2,Ii,il,i2) 
inc  =  (th2-thl)/Ii;  X  angle  increment 
theta  s  thl  +  inc*(0:li);  X  angles 
R  =  0; 

X  approximate  integral  by  Riemann  snm  at  each  angle;  take  care  of 
X  cases  when  denominator  conld  be  zero 
for  ii  3  i:li+l; 

if  il-i2*sin(theta(ii))  =*  0  ft  il+i2*sin(theta(ii))  ==0; 

R  =  R  +  2*(ob-oa)*inc; 
elseif  il-i2*sin(theta(ii))  ==  0; 

R  =R+((ob-oa)+(exp( j*ob*(il+i2*sin(theta(ii))))-  ... 
exp(j*oa*(il+i2*sin(theta(ii)))))/(j*(il+i2*sin(theta(ii)))))*inc; 
elseif  il+i2*sin(theta(ii))  ==  0; 
R=R+((exp(-j*oa*(il-i2*sin(theta(ii))))-  . . . 
exp(-j*ob*(il-i2*sin(theta(ii)))))/  . . . 
(j*(il-i2*sin(theta(ii))))+(ob-oa))*inc; 
else 

R  =R+((exp(-j*oa*(il-i2*sin(theta(ii))))-  ... 

exp(-j*ob*(il-i2*sin(theta(ii)))))/(j*(il-i2*sin(theta(ii))))+  . . . 
(exp(j*ob*(il+i2*sin(theta(ii))))-exp(j*oa3(il+i2*sin(theta(ii)))))/  . . . 
(j*(il+i2*sin(theta(ii)))))*inc; 
end 
end 

R  *  R/(4*pi); 

XRtdb.m 

X  thel  and  the2  are  the  lover  and  upper  angle  limits  of  the  source,  fl 

X  and  f2  are  the  lover  and  upper  frequency  limits  of  the  source,  and 

X  doc  is  the  interelement  spacing  divided  bt  the  speed  of  propagation 

function  [R]  =  Rtdb(K,L,thelIthe2lf 1, f2.doc) 

vb  =  2*pi*f2; 

va  =  2*pi*fl; 

thl  =  thel*pi/180; 

th2  =  th«2*pi/180; 

R  3  zeros (K*L, K*L) ; 

for  il  =  -L+1:L-1;  X  for  each  value  of  il,  compute  rov  and  column 
for  i2  3  0:K-1;  Xfind  rov  and  column  of  matrix  for  each  il 
rovd.^+l)  =  Rdb(va,vb,thl,th2,100,il,i2); 
col(l,i2+l)  =  Rdb(va,vb,thl,th2,100,il,-i2) ; 
end 

Rsub  =  toeplitz(col.rov) ; 
if  il  ==  0; 

R  =  R  +  kron(eye(L),Rsub);  X  add  matrices  to  main  diagonal 
elseif  il  >  0; 

X  add  matrices  to  diagonal  above  main  diagonal 
R  =  r  +  kron(diag(ones(l,L-il),il),Rsub); 
elseif  il  <  0; 

X  add  matrices  to  diagonal  belov  main  diagonal 
R  =  R  +  kron(diag(ones(l,L+il) ,il) ,Rsub) ; 
end 
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Xbfbbnfit.a 

ClMT 

K  *  16;  X  number  of  sensors 
L  =  5;  %  number  of  taps  par  sensor 

fl  =*  .2; 
f2  =  .4; 

doc  =  1;  %  interelement  spacing  divided  by  speed  of  propagation 

c_ang  a  -45;  X  constraint  angle  in  degrees 

sign_sq  a  .01;  X  noise  variance  or  poser 

sigl.sq  a  l;  %  variance  or  poser  of  signal  1 

sig2_sq  =  1;  X  variance  or  poser  of  signal  2 

R.thetal  =  mscmnes(K,L,45,f 1, f2.doc); 

R_theta2  a  mscmnes(K,L,-4S,f 1,12, doc); 

X  los  rank  approximation  of  R.thetal  and  R_theta2 

X 

[01, SI, VI]  =  svd(R_thetal) ; 
eigenl  a  svd(R_thetal) ; 

[Dl,per_enl]  a  getD(K,L, eigenl, .9999) 

U_L1  =  Ul(: ,1:D1); 

S_L1  =  S1(1:D1,1:D1) ; 

V_L1  =  Vl(: ,1:01); 

[U2,S2, V2]  =  svd(R_theta2) ; 
eigen2  =  svd(R_theta2) ; 

[D2,per_en2]  a  getD(K,L,eigen2, .9999) 

02  a  D1  X  ensure  matrix  sizes  are  same 
U_L2  =  U2(: ,1:02); 

S_L2  =  S2(1:D2,1:D2) ; 

V_L2  =  V2(: ,1:02); 

A1  a  o.Ll; 

Sigll  =  S_L1; 

A2  a  U_L2; 

Sig22  a  S_L2; 

X 

X 

X  another  say  to  get  square  roots  of  matrices 
XCvl.el]  =  eig (Sigll) ; 

XSigll_sr  =  vl*sqrt(el)*vl’ ; 

X[v2,e2]  =  eig(Sig22) ; 

XSig22.sr  =  v2*sqrt(e2)*v2»; 

XSigl2  =  Sigll_sr*Sig22_sr; 

X 

Sigl2  =  sqrt(Sigll)*sqrt(Sig22) ; 

Sig21  =  conj(Sigl2); 

XSigl2  =  zeros(size(Sigll)) ;  X  make  uncorrelated  for  testing  purposes 
XSig21  *  zeros (size(Sig22));  X  make  uncorrelated  for  testing  purposes 
R_ss  =  [A1  A2] * [Sigll  Sigl2;Sig21  Sig22]*[Al’ ;A2'3 ;  X  signal  covariance 
R_xx  =  R_ss  +  sign_sq*eye(K*L) ; 

XR_xx  =  sigl_sq*R_thetal  +  sig2_sq*R_theta2  +  sign_sq*eye(K*L) ;  X  for 
X  uncorrelated  testing 
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X 

X  coaqjute  eigenvector  conatrainta 
E.thetac  =  ascanev(K,L, c_ang, 11, f2.doc) ; 
eigen  *  avd(R_thetac) ; 

CD.per.en]  *  getO(K,L, eigen, .9998) 

[U , S , V]  *  avd(R_thetac); 

0_L  *  0(: .1:0); 

S_L  *  S(1:D,1:D); 

V.L  *  V(: ,1:D); 

C  *  V_L; 

li  =  nyfcen(K,L, c_ang.il, f2.doc); 
i  *  inv(S_L) *U_L ’  *ii ; 

*_o  *  C*inv(C**C)*f ; 

RIIV  =  inv(R.xx); 

*  =  RIIV*Oinv  ( C  * *RIIV*C) *i ; 
figured  ) 

plotvb2d(K,L, 1024, .2,*); 
hold 

plot «b2d(K,L, 1024, .225,*); 
plotvb2d(K,L, 1024, .25,w) ; 
plotwb2d(K, L, 1024. .275, ») ; 
plot *b2d(K,L, 1024, .3.*); 
plot wb2d(K,L, 1024, .325,*); 
plotvb2d(K,L, 1024, .35,*); 
plotvb2d(K ,L, 1024, .375,*); 
plot*b2d(K ,L , 1024 , .4 ,*) ; 
hold 

Xiigure(2) 

Xplot«b2d(K,L, 1024, .2,*); 

Xfigure(3) 

Xplot«b2d(K,L, 1024, .3,*); 

Xiigure(4) 

Xplot *b2d (K , L , 1024 , .4,*); 

XCvol.eol]  »  eig(Sl) ; 

XSl.ar  =  vol*aqrt(eol)*vol’; 

X[vo2,eo2]  =  eig(S2) ; 

XS2_sr  =  vo2*»qrt(eo2)*vo2’ ; 

Xior  ij  =  .2:. 025: .4; 

X  out_p  =  *'*Ul*diag(aqrt(Sl)) 

X  out_m  =  v'*02*diag(aqrt(S2)) 

X  out_p  =  *’*Ul*diag(aqrt(Sl))*l 

X  out_n  =  « ’ *U2*diag  C  sqrt ( S2 ) ) *rc*erp (- j  *2*pi*f  j  *del ) 
X  out_t  =  out_p  +  out_n 
X  aba (out _t) 

X  outdb_t  =  10*logl0(aba(out_t)) 

Xend 

Xout_p  »  *’*Ul*diag(aqrt(Sl))*l 

Xout_n  =  *'*02*diag(aqrt(S2))*rc*exp(-j*2*pi*.3*del) 

Xout_t  =  out_p  +  out_n 

Xaba(out_t) 

Xoutdb.t  =  10*logl0(aba(out_t)) 
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ont.pow  «  w'*E_aa*w 
aba(ont_pow) 

10*logl0(abs (out_pow) ) 
figured) 

Inaction  [1]  =  ■yicen(K,L, the, 11, f2.doc) 
wb  »  2*pi*f2; 
wa  *  2*pi*f 1; 
th  *  the*pi/180; 
z  =  doc*ain(th); 

E  *  zeroed,  K*L) ; 
ii  *  0; 

for  il  ■  (L-l)/2:-i:((L-l)/2)-L+i; 
for  12  =  (K-l)/2:-l:((K-l)/2)-K+l; 
ii  =  ii  +  1; 

fr(ii)  *  (l/(2*pi))*(«b*asinc(vb*(il+i2*x))-wa*aaiac(va*(il+i2*z))); 
and 
and 

f  =  ir.*; 

function  [answer]  =  nainc(z) 
top  =  length(x); 
for  ii  ■  l:top; 
if  x(ii)  ==  0; 

answer (ii)  =  1; 
also 

anavar(ii)  *  ain(x(ii))/x(ii); 
and 
and 

Xayfcanraf .n 

inaction  [f]  =  nyfcenref (K.L.tha, fl.f2.doc, P) 

*b  *  2*pi*f2; 
wa  *  2*pi*f 1; 
th  *  the*pi/180; 
x  =  doc*ain(th); 

E  =  zeros (1,K*L) ; 
ii  *  0; 

for  il  =  (L-i)/2:-i:((L-i)/2)-L+l; 
for  i2  =  (K-l)/2:-l:((K-l)/2)-K+l; 
ii  *  ii  +  1; 

fr(ii)  =  (l/(2*pi))*(wb*»inc(wb*(il+(P*sin(th))+i2*x))-  ... 
wa*nsinc(wa*(il+(P*sin(th))+i2*x))); 
and 
and 

f  »  fr.'; 

Xftdb.n 

Inaction  [f]  =  ftdb(K,L,thal, tha2.fl.f2.doc) 
wb  =  2*pi*f2; 
wa  =  2*pi*fl; 
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thl  »  th«l*pi/180; 
th2  ■  th«2+pi/180; 

E  *  z«roa(l,K*L); 

11  *  0; 

lor  il  *  (L-l)/2:-l:((L-l)/2)-L+l; 
lor  12  *  <I-l)/2:-l:«I-l)/2)-K+l; 
ii  ■  ii  ♦  1; 

lr(ii)  *  ldb(«a.nb,tbl.tb2, 100.il. 12)  ♦  idb(wa,»b,thl,tb2.100,il,-i2); 
•ad 
•ad 

1*1 r.»; 

Xldb.a 

Inaction  [1]  *  idb(oa,ob,thl ,th2, Ii,il, i2) 
inc  *  (th2-thl)/li; 
tb«ta  *  tbl  +  inc*(0:li); 

1*0; 

lor  ii  •  1:114-1; 
il  il+i2*ain(th*ta(ii))  **0; 

1*1+  2*(ob-oa)*iac; 

•ls« 

1  =  l+((axp(j*ob*(il+i2*ain(th«ta(ii))))-  ... 

•xp(j*oa*(il+i2*ain(th«ta(ii)))))/(j*(il+i2*ain(th«ta(ii))))+  . . . 
(•xp(-j*ob*(il+i2+8in(thata(ii))))-«xp(-j*oa*(il+i2*ain(th«ta(ii)))))/  . . . 
(-j*(il+i2*ain(th«ta(ii)))))*inc; 

•ad 

•ad 

I  *  l/(8+pi) ; 

%blbb4_23d.n 
clear 
L  *  6; 

K  *  16; 
doc  *  1; 
c_ang  *  18; 

II  *  .2; 

12  «  .4; 
aign.aq  =  .01; 
aigl.aq  *  10; 
aig2_aq  *  100; 

E_tb«tal  =  a>can«w(K,L, -17.8,11, 12.doc); 

R_th«ta2  =  ascan«n(K,L>-5.7S,ll,12,doc) ; 

El  *  E_intnr(K,L,100,-17.5,-6.7S,ll,12,«igl_aq,aig2_aq) ; 
nora(EI-(sigl_sq*E_th«tal  +  aig2_aq*R_thata2)) 

E_tb«tac  =  aacan«w(K,L,c_ang,ll,12,doc) ; 

•ig«n  *  ard(R_th«tac) 

Q),p«r_«n]  =  g«tO(X,L,*ig«n, .9999) 

XD  *c«il(((E-l)*(aia(c_aag+pi/180))+L-l)*2*(12-ll)+l) 

Xp«r_«n  =  ana(eig«n(l:D,l))/ann(«ig«n) 

E_zx  *  aigl_aq*E_tb«tal  +  aig2_aq*R_th«ta2  +  sign_sq*«y«(K*L) ; 

E_zxa  *  El  +  aiga_aq*«y«(X*L); 
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[U , S , V]  =  «Yd(R_thetac) ; 

0.L  *  U(: ,1:D); 

S_L  *  S(1:D,1:D) ; 

¥_L  *  V( : ,1:D) ; 

C  =  V_L; 

f  i  =  ayfcen(K,L,e_aiig, fl.f2.doc); 
f  *  imr(SJL)*U_L’efi; 
w_o  *  C*ixnr(C’*C)*f ; 

CC  *  [C  w_o) ; 

DJU.SS.W]  »  svd(CC) ; 

C_n  =  UU(:  ,D+1:K*L); 

C_n'*C 

C_n’*w_o 

D 

per.en 

w_n  =  inv (C_n ’ *R_xx*C_n) *C_n ’ *R_xx*w_ 
w2  *  w_o  -  C_n*w_n; 

RIIV  a  inv(R_xx) ; 

RUVn  =  inv(R_xxn) ; 
wl  a  RHV*C*inv(C'*RIIV*C)*f ; 
mi  =  RI«Vn*C* inw  C  C  *  *RIIVn*C) *f ; 
Xfigured) 

Xplotfftwb2(K,L,512, .3,wl,wn) ; 
Xfigure(2) 

Xplotfftwb2(K,L,512. .3,wl,w2); 

Xf igure(3) 

Xplot  vb2d2 (K , L , 1024 , .3,wl,wn) ; 
Xfigure(2) 

Xplot wb2d2 (K , L , 1024 , .3,wl,w2); 
figure(l) 

Xplotwb2d(K,L,1024, .3,w_o); 
plotwb2d(K,L,1024, .2,wl) ; 
hold 

plotwb2d (K ,  L , 1024 , .225,wl); 
plotwb2d(K,L, 1024, .25,wi) ; 
plotwb2d(K,L, 1024, .275,wl); 
plotwb2d (K , L , 1024 , . 3 , wl ) ; 
plotwb2d(K,L, 1024, .326,wl); 
plotwb2d(K , L , 1024 , .35,wl) ; 
plotwb2d(K,L, 1024,. 375, wl); 
plot vb2d(K,L, 1024, .4,wl) ; 
hold 

f igure(2) 

plotwb3d(K,L,32,wl) ; 

Xfigure(2) 

Xplotvb2d2(K ,L, 1024, .3,ul ,w2) ; 

Xf igure(3) 

Xplotgf 1 (K , L , 1024 , c_ang , w_o) ; 
figure (4) 

plotgf 1(K,L, 1024, -17.5, wl) ; 
hold 


plotgll(K,L, 1024,-5. 75, wl); 
hold 

ligora(S) 

plotpl 1 (K , L , 1024 , c.ang , w_o) ; 
ligora(l) 

Xplot«b3d.a 

function  [B3  =  plot«b3d(K,L,I,«) 
ouga  *  2*pi/I*(0:I)  -  pi; 
thata  =  pi/I* (0:1)  -  pi/2; 

E  *  zaros(I+l,I+l) ; 
d  *  zaros(K*L+(I+l) ,1+1) ; 
lor  il  =  1:1+1; 
lor  i2  =  1:1+1; 

d(l+(il-l)*K*L:il*K*L,i2)=kron(axp(-j*oaaga(il)) .‘([0:L-l3 .  *) 
,axp(j*oaaga(il)*sin(thata(i2))) .“([(K-l)/2:-l: ((I-l)/2)-I+l3 . *)) ; 
and 

R(il, : )  =  *’*d(l+(il-l)*K*L:il*K*L, :); 
and 

aashc (tbata* 180/pi ,oaaga, 10*logl0 (aba (E) . *2)) 
zlabalC’DOA  in  dagraaa' ) 
ylabal( 'formalized  1  in  rad') 
zlabal('Gain  in  dB*) 

Xblbbl4_9.a 
claar 
L  *  8; 

K  =  15; 

11  =  .3; 

12  *  .4; 

dOC  =  1; 

c_ang  *  6.74; 

R_thatac  =  B8can*v(K,L,c_ang,ll,12,doc) ; 

E.thetal  =  aacanaw(K,L,23. 58, 11,12, doc) ; 

E_thata2  =  ascanav(K,L,36.87,ll,12,doc); 
aigan  =  avd(R_thatac) 

[D,par_an3  =  gatD(K,L, aigan, .9999) 
aign_sq  =  .01; 
aigl.aq  *  10; 
aig2_aq  =  100; 

R_xx  =  sigl_sq*R_thetal  +  sig2_sq*E_thata2  +  sign_aq*aya(K*L) ; 
CO.S.y)  =  avd(E_thatac) ; 

U_L  =  0(:,1:D); 

V_L  =  V(: ,1:0) ; 

S_L  =  S(1:D,1:D); 

C  *  V_L; 

li  =  aylcan(K,L,c_ang,ll,12,doc); 

1  =  inv (S_L) *0_L ’ *li ; 
w_o  =  C*inv(C**C)*l; 

El  =  inv(R.xx); 

»  *  EI*C*inr(C'*EI*C)*l; 
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2igur«(l) 

plot*b2d(K,L, 1024. . 3 ,«) ; 
hold 

plot«b2d(E,L, 1024. .326, ») ; 
plot«b2d(K , L , 1024 , . 35 , ») ; 
plot«b2d(K . L , 1024 . . 375 , ■) ; 
plotvb2d(K,L, 1024, .4,*); 
hold 

figured) 

plotgf 1 (K , L , 1024 , c.ang , w.o ) ; 
figured) 

plotpf  1  (K ,  L ,  1024 ,  c.ang ,  w.o) ; 
figure (4) 

plotwb2d(K,L, 1024, .3, w.o)  ; 
hold 

plotwb2d(K,L,1024, .325,u_o); 
plotvb2d(K,L, 1024, .35, w.o); 
plotwb2d(K,L, 1024, .375, w.o); 
plotwb2d(K ,L , 1024, .4,»_o) ; 
hold 

figure(l) 

%bf bbadddb . m 
clear 
L  =  5; 

K  *  16; 
cod  =  1; 
c.ang  =  0; 
c.angl  =  -6; 
c_ang2  =  0; 
fl  =  .2; 
f2  =  .4; 
sign.sq  =  .01; 
aigl.sq  =  1; 
sig2_sq  *  1; 

R.thetal  =  Rtdb(K,L,30,35,fl,f2,l) ; 

%R_theta2  =  Rtdb(K,L,40.45,fl,f2,l); 

R.thetac  =  Rtdb(K, L, c.angl, c_ang2,fl,f2,l) ; 

XR.thetac  =  aacamew(K, L, c.ang, f 1,12,1) ; 
eigen  =  svd(R.thetac) 

X[D,per_en]  =  getD(K,L, eigen, .87) 

D  =ceil(((K-l)*(sin(c_ang2*pi/180)-ein(c_angl*pi/180))+L-l)*2*d2-fl)+l) 
per.en  =  aum( eigen (1 :D, 1)) /sun (eigen) 

%R_xx  =  aigl.aqeR.thetal  +  aig2_aq*R_theta2  +  sign_aq*eye(K*L) ; 

XR.zxn  =  RI  +  aign_aq*eye(K*L) ; 

R. xz  =  aigl_aq*R_thatal  +  aign_aq*eye(K*L) ; 

XR.xz  ■  sign_aq*eye(K*L) ; 

[U , S , V]  =  avd(R.thetac) ; 

0.L  =  U(: ,1:D); 

S. L  =  S(1:D,1:D); 

V.L  =  V(:,1:D); 
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c  *  V_L; 

fi  *  ftdb(K,L,c_angl,c_ang2,fl,12,l); 

Xfi  =  ayfcen(K,L,c_ang,fl,12,l); 

I  »  in* (S_L) *U_L ’  *1  i ; 

»_o  *  C*in*(C’*C)*f ; 

El IV  *  in»(R_xx) ; 

XRIIVn  «  in*(R_xxn) ; 

v  »  RHV*C*inT  (C ' *RI*V*C ) *1 ; 

X»n  *  RlIVneOin*  (  C  *  *RIIVn*C)  *1 ; 

Xfigured) 

Xplotf 1 tvb (K , L , 6 12 , .3.*); 

Xfigure(2) 

Xplotvb2d(K,L, 1024, .2,*); 

Xhold 

Xplot*b2d(K , L , 1024 . . 225 , ») ; 

Xplotvb2d(K , L , 1024 , .26,*) ; 

Xplot*b2d(K ,L, 1024 , .275,*) ; 

Xplotvb2d(K,L, 1024, .3,*); 

Xplot«b2d(K , L , 1024 ,. 326 ,«) ; 

Xplot«b2d(K , L , 1024 , .36,*) ; 

Xplotvb2d (K , L , 1024 , .376,*); 

Xplotvb2d(K,L, 1024, .4,*); 

Xhold 

figured) 

plots! 1 (K , L , 1024 , ( ( c.ang l+c_ang2) /2 ) , *) ; 
figured) 

plotpf 1(K,L, 1024, ((c_angl+c_ang2)/2) ,«) ; 
ligur*(l) 

Xbfbbrefnn.B 
clear 
K  =  16; 

L  =  16; 

II  =  .2; 

12  =  .3; 
cod  *  1; 
c_ang  =  -46; 

P  -  9; 

rc  =  .6; 
sign_sq  =  .01; 

R_ll  =  ascanev(K,L,46,fl,12,cod); 

R_21  =  asca21h(K,L,46,fl,f2,cod,P); 

R_12  =  BScal2h(K,L,46,fl,12,cod,P); 

R_22  =  bbcb22(K,L, 45, 11,12, cod); 

XR.22  =  ■scane«(K,L, -46, 11,12, cod); 

R_xx  =  R_ll  +  rc*R_21  +  rc*R_12  +  (rc.‘2)*R_22  +  sign_8qe*ye(K*L) ; 
R_thetac  >  ascBnev(K,L,c_asg,ll,12,cod) ; 
eigen  -  »vd(R_thetac) ; 

[D.per.en]  =  getD(K,L, eigen, .9999) 

[0,S,V]  =  svd(R.thetac) ; 

U_L  =  U(:,1:D); 
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S_L  *  S(1:D,1:D); 

V_L  *  V(: ,1:D) ; 

C  *  V_L; 

fi  *  ayfcearef (K,L,c_aag,f 1,12, cod, P) ; 
f  *  ia»(S_L)*U_L'*fi; 
w_o  *  Oiav(C»*C)*f ; 

RI1V  »  inv(R_xx); 

v  =  RI  W*C*iax  (C '  *RIIV*C  )  *t ; 

ligar«(l) 

plotwb2d(K,L, 1024,  .2,e) ; 
hold 

plot vb2d(K,L, 1024, .225,u) ; 
plotwb2d(K,L, 1024, .25,«) ; 
plotvb2d(X , L , 1024 . . 276 , ■) ; 
plotvb2d(K , L , 1024 , .3,w); 
plot«b2d(K,L, 1024, .326,*) ; 
plotwb2d(K,L, 1024, .36,*); 
plotvb2d(K , L , 1024 , .376,*) ; 
plot*b2d(K , L , 1024 , .4,u); 
hold 

%plotvb2d(K,L, 1024, ,3,»); 
out_p  =  w'*R_xx*u 
10*logl0(out_p) 

Xoutputpounbq . ■ 
clear 
d  =  1; 
n  *  16; 
lambda  =  2; 
sigmaa.sq  =  le-6; 
rc  ■  .26; 

P  =  8; 

del  =  2*P* . 707 ; 

R_aa  =  [1  rc*exp( j*2*pi/lambda*del*d) ;  ... 

rc*exp(-j*2*pi/lambda*del*d)  rc. ‘2] ; 
did  =  exp(j*2*pi*d/laabda*(0.707)).‘([(a-l)/2:-l:((a-l)/2)-a+l]  .  ’); 
dii  =  exp(j*2*pi*d/lambda*(-0.707)) .*([(a-l)/2:-l: ((a-l)/2)-a+l] . ’) 
ds  =  exp(j*2*pi*d/laabda*(-0.707)) .  “([(a-l)/2:-l : ((a-l)/2)-a+l]  . ') ; 
c  *  da; 
f  =  l; 

u_o  =  c*iav(c'*c)*f ; 
cc  =  [c  w_o] ; 

[U.S.V]  =  svd(cc) ; 

C_a  =  U(:,2:a); 

C_a’*c 

C_a ’ *u_o 

H  =  [did  dii] ; 

R_xx  =  H*R_aa*H '  +  sigmaa_sq*eye(a) ; 

RI  =  iav(R_xx) ; 

wl  =  RI*c*iav(c’*RI*c)*f ; 

w_n  =  iav(C_a’*R_xx*C_a)*C_a'*R_xx*w_o; 
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w2  =  w_o  -  C_n*w_n; 
op  =  zeros(181,l); 
opdb  =  zeros(l81 , 1) ; 
lor  ang  *  -90:0; 
angr  =  ang*pi/180; 

di  =  exp(-j*2*pi*d/lambda*sin(angr))  .*  ([0:n-l] . ') ; 

H  *  [di] ; 

Ex  =  H*H'  +  sigman_8q*eye(n) ; 
op(ang+91)  *  w_o’eRx*w_o; 
opdb(ang+91)  =  10*logl0(op(ang+91)) ; 
end 

lor  ang  =  1:90; 
angr  =  ang*pi/180; 
rc  »  (angr‘2)/((pi/i)"2) ; 
del  =  2*P*sin(angr) ; 

R_aa  =  [1  rc*exp( j*2*pi/lambda*del*d) ;  ... 

rc*exp(-j*2*pi/lambda*del*d)  rc.*2] ; 
did  =  exp(j*2*pi*d/la»bda*sin(angr)).*([(n-l)/2:-l:((n-l)/2)-n+l]  .  *); 
dii  =  exp( j*2*pi*d/lambda*sin(-angr) ) .*([(n-l)/2:-l: C(n-i)/2)-n+l] . '); 
H  =  [did  dii] ; 

Ex  =  H*R_aa*H’  +  8igman_sq*eye(n) ; 
op(ang+91)  =  w_o’*Ex*w_o; 
opdb(ang+91)  =  10*logl0(op(ang+91)) ; 
end 

angles  =  [-90 : 90] ; 

plot (angles , 10*logl0(abs(op) ) ) 

xlabel('DOA') 

ylabel( ’Output  Power  in  dB  with  Quiescent  Heights') 

Xoutputpownbw . m 
clear 
d  =  1; 
n  =  16; 
lambda  =  2; 
sigman_sq  =  le-6; 
rc  =  .25; 

P  =  8; 

del  =  2*P*.707; 

E_aa  =  [1  rc*exp( j*2*pi/lambda*del*d) ;  ... 

rc*exp(-j*2*pi/lambda*del*d)  rc.‘2]; 
did  =  exp (j*2*pi*d/lambda* (0.707)) . *( [(n-l)/2:-l : ((n-l)/2)-n+l]  . ') ; 
dii  =  exp ( j*2*pi*d/lambda* (-0.707)) . *([(n-l)/2:-l : ((n-l)/2)-n+l] . ’) ; 
ds  =  exp(j*2*pi*d/lambda* (-0.707)) . *([(n-l)/2:-l: ((n-l)/2)-n+l] . ') ; 
c  =  ds; 

1  =  1; 

w_o  =  c*inv(c'*c)*l; 
cc  =  [c  w_o] ; 

[U , S , V]  =  svd(cc) ; 

C_n  =  U( : ,2:n) ; 

C_n’*c 

C_n’*w_o 
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H  *  [did  dii] ; 

R_xx  =  H*R_aa*H’  +  sigman_8q*eye(n) ; 

RI  *  inx(R_xx); 

wl  =  RI*c*inv(c’*RI*c)*l; 

w_n  =  inv(C_n ’ *R_xx*C_n) *C_n ’ *R_xx*b_o ; 

*2  =  w_o  -  C_n*w_n; 
op  =  z«ros(181,l); 
opdb  =  z«ros(l81 , 1) ; 
lor  ang  =  -90:0; 
angr  =  ang*pi/180; 

di  *  exp(-j*2*pi*d/laabda*sin(angr))  ([0:n-l].'); 

H  =  [di]; 

Rx  =  H*H’  +  8igman_sq*eye(n) ; 
op(ang+91)  =  w2'*Rx*w2; 
opdb(ang+91)  =  10*logl0(op(ang+91)); 
and 

lor  ang  =  1:90; 
angr  =  ang*pi/180; 
rc  =  (angr“2)/((pi/2)“2); 
dal  =  2*P*sin(angr) ; 

R_aa  a  [1  rc*exp( j*2*pi/laabda*dal*d) ;  ... 

rc*exp(-j*2*pi/laabda*del*d)  rc.”2] ; 
did  =  exp(j*2*pi*d/laiabda*ain(angr)).‘([(n-l)/2:-l:((n-l)/2)-n+l]  .  ’); 
dii  =  exp(j*2*pi*d/lanbda*sin(-angr)).*([(n-l)/2:-l:((n-l)/2)-n+l]  . '); 

H  =  [did  dii] ; 

Rx  =  H*R_aa*H ’  +  signan_aq*eye(n) ; 
op(ang+91)  =  v2  1  *Rx*w2 ; 
opdb(ang+91)  =  10*logl0(op(ang+91)) ; 
and 

anglas  =  [-90 : 90] ; 
plotCangles, 10*logl0(abs(op))) 
xlabel(’DOA’) 

y label ( ’Output  Power  in  dB  with  Adapted  Weights ’ ) 

%out putpowq . m 
K  a  16; 

L  a  25; 

11  =  .3; 

12  =  .4; 
doc  =  1; 
c_ang  =  -45; 
ang  =  46; 

P  =  9; 
rc  =  0.25; 
sign.sq  =  .01; 

R_ll  =  nscmnew(K,L,ang,ll,12,doc); 

R_21  =  B8CB21h(K,L,ang,ll,12,doc,P); 

R_12  =  mscail2h(K,L,ang,ll,12,doc,P); 

R_22  a  mscm22(K,L,ang,ll,12,doc); 

R.thetac  =  mscmnew(K,L,c_ang,f 1 ,12,doc) ; 

R_xx  =  R_ll  +  rc*R_21  +  rc*R_12  +  (rc.“2)*R_22  +  sign_sq*eye(K*L); 


eigen  =  svd(R.thetac) ; 

[D,per_en]  =  getD(K,L, eigen, .9999) ; 

[U,S,V]  =  svd(R_thetac) ; 

U_L  =  U(: ,1:D) ; 

S_L  =  S(1:D,1:D); 

V_L  *  V(:,1:D); 

C  =  V_L; 

li  =  ayl cenr el (K, L , c_ang, f 1, f2.doc, P) ; 

I  =  inv (S_L) *U_L 1 *1 i ; 

»_o  =  C*inv(C’eC)*l; 

CC  =  [C  w_o] ; 

[UU.SS.W]  =  svd(CC) ; 

C_n  =  00( : ,D+1 :K*L) ; 

w_n  =  inv (C_n ’ *R_xx*C_n) *C_n ’ *R_xx*w_o ; 

*2  =  »_o  -  C_n*w_n; 

RIIV  =  inv(R.xx); 

«  =  RIIV*Oinv(C ' *RIIV*C) *1 ; 
op  =  zeros (181,1) ; 
opdb  =  zerosC 181,1) ; 
lor  ang  =  -90:0; 

R  =  ascanev(16,L,ang,ll,12,doc); 

Rx  =  R  +  sign_sq*eye(K*L) ; 
op(ang+91)  =  s_o’*Rx*w_o; 
opdb(ang+91)  =  10*logl0(op(ang+91)) ; 
end 

lor  ang  =  1:90; 
angr  =  ang*pi/180; 
rc  *  (angr"2)/( (pi/2) *2) 

R_ll  =  ascnnev(K,L,ang,ll, 12.doc); 

R_21  =  B8ca21h(K,L,ang,ll,12,doc,P); 

R_12  =  ascal2h(K,L,ang,ll ,12, doc, P); 

R_22  =  B8ca22(K,L,ang,ll,12,doc); 

Rx  =  R_ll  +  rc*R_21  +  rc*R_12  +  (rc“2)*R_22  +  sign_sq*eye(K*L) ; 
op(ang+91)  =  w_o’*Rx*w_o; 
opdb(ang+91)  =  10*logl0(op(ang+91)) ; 
end 

opaax  =  aax(op) ; 

opnorn  =  op/ opaax, 

angles  =  [-90 : 90] ; 

plot (angles , 10*logl0( opaax) ) 

xlabel(’DOA’) 

ylabel( 'formalized  Output  Power  in  dB  with  Quiescent  Weights’) 

%outputpo««.a 
K  =  16; 

L  =  25; 

II  *  .3; 

12  =  .4; 
doc  *  1; 
c_ang  =  -45; 
ang  =  45; 
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p  =  o; 

rc  =  0.25; 
sign.sq  =  .01; 

R_ll  =  nscnneeCK.L, ang, 11, 12.doc)  ; 

R_21  =  nsca21h(K,L>ang,ll,12,doc,P); 

R_12  =  nscal2MK,L, ang,  11, 12.doc, P) ; 

R_22  =  ascn22(K,L,ang,ll ,12, doc) ; 

R_thetac  =  a8camev(K,L,c_ang, 11,12, doc) ; 

R_xx  =  R_ll  +  rc*R_21  +  rc*R_12  +  (rc.‘2)*R_22  +  sign_sq*eye(K*L) 
eigen  =  evd(R_thetac) ; 

[D,per_en]  =  getD (K.L, eigen, .9999) ; 

[U,S,V]  =  svd(R_thetac) ; 

0_L  =  U(:,1:D); 

S_L  =  S(1:D,1:D); 

V_L  =  V ( : ,1:D); 

C  =  V_L; 

li  =  mylcenrel(K,L,c_ang,ll,12,doc,P) ; 

1  =  inv(S_L) *U_L ’  *li ; 
w_o  *  C*inv(C’*C)*l ; 

CC  =  EC  w_o] ; 

[UU.SS.W]  »  svd(CC) ; 

C_n  =  UU(: ,D+1:K*L); 

w_n  =  inv (C_n ’ *R_xx*C_n) *C_n  ’ *R_xx*w_o ; 

*2  =  w_o  -  C_n*®_n; 

RIIV  s  inv(R_xx) ; 

w  =  RI*V*C*inv (C ’ *RIIV*C)*1 ; 

op  »  zeros (181,1) ; 

opdb  =  zeros(181,l); 

lor  ang  -  -90:0; 

R  =  nscnnew(16,L,ang,ll,12,doc); 

Rx  =  R  +  sign_sq*eye(K*L) ; 
op(ang+91)  =  w’*Rx*w; 
opdb(ang+91)  =  10*logl0(op(ang+91)) ; 
end 

lor  ang  =  1:90; 
angr  =  ang*pi/180 ; 
rc  =  (angr“2)/((pi/2)“2) ; 

R_ll  =  ascnnev(K,L,ang,ll,12,doc); 

R_21  =  asca21h(K,L, ang, 11, 12.doc ,P); 

R_12  =  nscnl2h(K,L,ang,ll,12,doc,P); 

R_22  =  B8cn22(K,L,ang,ll,12,doc) ; 

Rx  *  R_ll  +  rc*R_21  +  rc*R_12  +  (rc‘2)*R_22  +  sign_sq*eye(K*L) ; 
op(ang+91)  =  b’*Rx*b; 
opdb(ang+91)  =  10*logl0(op(ang+91) ) ; 
end 

opaax  =  Bax (op) ; 

opnon  =  op/opaax, 

angles  =  [-90 : 90] ; 

plot (angles, 10*logl0(opaax)) 

xlabel(’DOA’) 

y label ( ’lormalized  Output  Power  in  dB  with  Adapted  Weights’) 
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